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ABSTRACT 


The  phase-plane  method  is  a  graphical  method  for 
linear  and  non-linear  second-order  systems.   There  are 
many  techniques  for  obtaining  the  phase  portrait  which 
consists  of  a  number  of  phase  trajectories  on  the  phase- 
plane.   From  the  summary  and  discussion  the  isocline 
method  is  a  most  general  and  useful  method.   The  only 
difficulty  is  in  the  labor  required.   Solution  by  digi- 
tal computer  reduces  this  difficulty  and  makes  the  iso- 
cline method  a  much  more  useful  method  in  non-linear 
systems  analysis  and  design  applications. 
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I.   INTRODUCTION 

A.  CONDITIONS  FOR  LINEARITY  AND  DEFINITION 
OF  NON-LINEARITY 

It  is  useful  to  define  linearity  and  then  any  system, 
component,  or  phenomenon,  which  does  not  satisfy  this  de- 
finition will  be  considered  non-linear.   The  simple  mean- 
ing of  linearity  is  that  if  the  relationship  between  output 
y  and  input  x  is  plotted  the  result  is  a  straight  line.   Or 
that  if  the  output  is  proportional  to  the  input  the  process 
must  be  linear.   So  we  can  simply  say  the  essential  condi- 
tion for  linearity  is  that  the  law  of  superposition  applies, 
and  linearity  is  a  specification  which  therefore  limits  the 
field  of  activity.   Non-linearity  is  a  "non-specification" 
and  the  field  is  unbounded,  linearity  being  a  particular 
part  of  it. 

B.  CLASSES  OF  NON-LINEARITY 

Non-linearity  in  physical  systems  cannot  be  avoided, 
and  in  fact  can  sometimes  be  used  to  improve  the  system's 
performance.   The  non-linearities  encountered  in  control 
systems  may  be  classified  as  accidental  non-linearities 
and  intentional  non-linearities.   The  accidental  class  con- 
sists of  those  non-linearities  which  are  inherent  in  the 
components  used  in  the  system.   Some  of  these  are  coulomb 
friction,  stiction,  backlash,  binding,  dead  zone,  satura- 
tion, etc.   Their  existence  is  accidental  in  the  sense  that 


15 


they  are  not  deliberately  inserted  in  the  system.  Inten- 
tional non-linearities  are  those  which  are  introduced  de- 
liberately, either  to  obtain  desired  performance  features 
or  for  weight,  space,  or  economic  advantages. 

Clearly  the  simplest  class  of  non-linearity  is  the  in- 
stantaneous, a  member  of  which  only  changes  its  input  u(t) 
to  some  output  x(t)  without  any  delay,  the  ratio  x(t)/u(t) 
being  a  function  of  u(t).   Thus  x(t)  =  f[u(t)]  and  it  can 
be  represented  in  schematic  form  (see  Fig. (1-1)),  having  a 
transfer  function  f(z)/z  which  is  independent  of  time. 
Saturating  amplifiers  and  relays  are  examples  of  instan- 
taneous non-linearities.   In  practice,  of  course,  such  a 
non-linearity  will  almost  always  be  associated  with  some 
linear  (or  non-linear)  weighting  function  resulting  from 
inertia  and  damping  as  in  Fig.  (1-2), so  that 


x(t)  =  f[u(t)  ] 


(1-D 


In  general,  the  instantaneous  non-linear  characteristic  is 
combined  with  inertia  and  damping  in  the  same  component,  in 
which  case  it  cannot  be  separated  as  in  Fig. (1-2) . 


u(t) 


x(t) 


f  (z) 
Fig.  1-1.   A  transfer  function    jj  '  which 

is  independent  of  time. 
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u(t) 

^-- 


y(t) 


Fig.  1-2.   The  instantaneous  non-linear 

characteristic  is  combined  with 
inertia  and  damping  cannot  be 
separated  as  in  figure. 

By  convolution  integral 


y(t)  =  j  x(T)h(t-T)dT  ,  and 


changing  the  variable  gives 


y(t)  = 


x(t-T1)h(T1)dr1 


Assume 


x  <  0    for    t  <  0 


so  that 


x(t-T1)  =  0    for    T1  ^  t 
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Then 


y(t)  =  J  x(t-T1)h(T1)dT. 
o 


l'"w.l'u,l 


f[u(t-T1) ]h(T1)dT1 


r»t 


KCuCt-T^^)  ,  T1'\dT1 


(1-2) 


K[u ( t-T, ) , T, ]  is  called  the  Kernel  function  and  is  seen  to 
depend  on  both  the  dynamics  of  the  system,  f(z)/z  and  h(t), 
and  the  characteristics  of  the  input  process,  u(t). 

Zadeh  has  classified  non-linear  systems  on  the  basis 
of  multiway  integrals  and  Kernel  functions.   Thus  those 
that  can  be  represented  by  a  single-way  integral  belong  to 
the  first  class,  r?,  ,  which  therefore  includes  the  subclass 
of  all  linear  systems  and  also  the  degenerate  zero-order 
class  of  instantaneous  non-linearities. 

Then  the  first-order  non-linear  differential  equation 
is 

dx 


H  +  Fq(x)  =  P1[u(t)] 


or 


If  =Fl[u(t)] 


F  (X) 


Integrating  gives 


x(t)  =    K[u(t-T) , T]dT 
o 
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Hence,  the  second-order  non-linear  differential  equation  is 


f  (x,  x,  x)  =  u(t) . 

Integrating  twice  gives 

x(t)  = 


K 


o  o 


[uCt-T^,  U(t-T2) ,T1#T2]dT1dT2 


So,  the  second  class,  r\~,    contains  all  those  for  which 
the  input/output  relationship  can  be  expressed  as  a  2-way 
integral  and  this  includes  the  degenerate  cases  belonging 
to  class  r\   . 

The  nth  class,  T?  ,  contains  all  those  systems  for 
which  the  input/output  relationship  is  a  n-way  integral,, 
which  means  that  they  can  be  represented  by  nth  order  non- 
linear differential  equations. 

Zadeh' s  classification  is  an  important  property  in 
considering  random  processes  in  the  presence  of  noise. 

C.    NON-LINEAR  ANALYSIS 

Many  methods  used  in  the  analysis  of  non-linear  prob- 
lems are  effective  only  for  limited  types  of  problems.   A 
few  of  them  have  been  developed  which  have  more  general 
applicability  to  non-linear  control  systems: 
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1.  Linear  approximations. 

2.  Graphical  method  -  The  phase-plane. 

3.  Piecewise  linear  approximation. 

4.  Numerical  methods. 

5.  Describing  function. 

6.  The  analog/digital  computer. 

D.    PHASE-PLANE  ANALYSIS  OF  NONLINEAR  SYSTEMS 

The  phase-plane  method  is  a  graphical  method.   There 
are  many  techniques  for  obtaining  the  phase  portrait  which 
consists  of  a  number  of  phase  trajectories  on  the  x  vs  x 
plane  or  in  some  plane  related  to  this  plane  through  a  con- 
venient transformation  of  variables.   After  the  phase  por- 
trait is  obtained,  by  means  of  it,  analysis  of  the  follow- 
ing aspects  and  items  is  possible: 

1.  Computation  of  the  exact  transient  response  for  a 
given  set  of  initial  conditions. 

2.  Prediction  of  the  type  of  transient  response 
(oscillatory  or  monotonic)  with  a  certain  range  of  initial 

conditions. 

3.  Estimates  of  the  amount  of  overshoot  in  response 
to  a  step  input  of  given  magnitude. 

4.  Prediction  of  the  system  stability. 

5.  Prediction  of  the  existence  and/or  amplitude  of 
limit  cycles. 

6.  Prediction  of  the  existence  of  various  modes  of 
operation. 
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7.  Study  of  the  effect  of  a  given  type  of  non- 
linearity  on  the  performance  of  the  system  as  compared 
with  a  linear  system,  with  a  system  having  a  different 
type  of  non-linearity,  or  with  a  system  which  is  non-linear 
in  the  same  sense  but  in  different  degree. 

8.  Attempts  to  answer  questions  such  as: 

a.  Does  this  non-linearity  improve  performance? 

b.  What  changes  can  be  made  to  suppress  the  limit 
cycle? 

c.  How  can  we  adjust  this  system  to  make  it  meet 
specifications? 

As  stated  above,  we  can  understand  a  graphical  method; 
it  is  usually  relatively  simple  to  utilize,  and  may  be 
especially  useful  as  an  exploratory  tool  when  a  non-linear 
equation  is  first  being  attacked.   There  are  many  ways  to 
obtain  phase  trajectories  on  the  phase-plane.   In  the  next 
chapter,  the  most  useful  methods  will  be  summarized  and 
compared.   The  advantages  and  disadvantages  of  the  various 
methods  will  be  pointed  out. 
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II.   SUMMARY  OF  USEFUL  GRAPHICAL  METHODS 
FOR  CONSTRUCTING  PHASE  TRAJECTORIES 
ON  THE  PHASE-PLANE 

A.    ISOCLINE  METHOD 

1.   Method-description 

The  free  motion  of  any  second-order  non-linear 
system  can  be  described  by  an  equation  of  the  form: 

x  +  $(x,x)x  +  >£(x,x)x  =  0 

where 

dx 

dt 

d  x 

x  =  — 2 

dt 


Let 


then 


y  -  f  (2-1) 


*  =  dt  =  ~  *<*'¥>¥  "  ^X'Y)X  <2~2) 


Equation  (2-2)  divided  by  (2-1)  yields 


dy_       <E(x,y)y  +  >£(x,y)x 
dx  y 


Let 


-^-  =   M   where  M  is  a  constant, 
dx 
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Then 

My  =  $(x,y)y  +  >£(x,y)x  (2-3) 

Equation  (2-3)  is  the  equation  of  isoclines. 
NOTES:    a.   If  $(x,y)  and  ^(x,y)  are  constants,  then 
isoclines  are  straight  lines.   Otherwise,  isoclines  will 
be  curves. 

b.  For  a  first-order  equation: 

dt  =  f(x'fc) 
Isoclines  will  be   f(x,t)  =  M. 

c.  Equations  more  than  second-order  need  some 
transformation  technique,  so  this  method  is  good  only  for 
second-order  systems. 

d.  After  obtaining  the  isoclines,  phase  trajec- 
tories are  plotted  on  the  phase-plane  as  follows: 

M, ,  Mp,  M_  represent  the  slope  of  each 

isocline.   Start  from  same  initial  condition  point,  as  in 
Fig.  2-1  point  A.   From  point  A  draw  dotted  lines  parallel 
to  M..  ,  M_.   These  dotted  lines  intersect  the  M~  isocline 
at  a  and  be.   B  is  middle  point  of  ab.   Continuing  in  this 
way,  obtain  C,  D  ...  points.   Connect  A,  B  ...  to  obtain 
a  phase  trajectory  from  the  given  initial  condition  A  on 
the  phase-plane. 
2.   Example 

Consider  a  viscous-damped  servo 

d29       d9 

J  —t~-  +   F  -r-£  +  K0   =0. 
dt       dt       o 
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-»-x 


Fig.  2-1.   Construction  of  a 
trajectory  using 
isoclines. 
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Dividing  through  by  J: 

d2e     _   d0 

— =S  +  £  —P-  +  |  0   =0.  (2-4) 

,,2     J   dt    J   o 
dt 


Substituting 


"k  =  J  J 


and 

F 


c   = 


2^/KJ 
into  equation  (2-4) 


d20  d0  9 

r2-  +  2co)     -T7T-  +    U>      6     =    0. 

,,2  n   dt  no 

dt 


Let 


dt 
then 


°   =    z 


-,  ,      ^  d0     x        ,      d  e 

dz  d     /         o    \        1  o 


d0  d0      V   dt      J  "   z       ,.2 

o  o  dt 


so 

dz 


z    = 


+    2CU>   Z    +     05       0        =    0 


d0      T   ^"h"    T    "h      O 


Let 

dz 

M  =   dfl- 
o 
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Isocline  equations  will  be 

ZM  =  -2cb)   z  -  oj  20  .  (2-5) 

n     n   o  \      ->  i 

If  c  =  0.25  (linear  servo  system) 

and  0)  =1, 
n 

equation  (2-5)  will  be  ZM  =  -0.5z  -0 
^  o 

If  M  is  equal  to  2,1,0,  -0.5,  -1,  -2  we  obtain  straight 
radial  lines  (isoclines)  as  in  Fig.  2-2. 

Assume  an  initial  condition  represented  by  A  in  the 

figure,  namely  an  output  position  0   =  -2.7  and  a  position 

dBo  ° 

rate  of  change  -=77-  =  0,  we  can  plot  the  phase  trajectory  as 

in  Fig.  2-2. 

3 .   Discussion 

a.  From  Fig.  2-2  it  is  seen  that  by  plotting  some 
isoclines  on  the  phase-plane  (more  or  less  depending  on  de- 
sired accuracy)  then  any  phase  trajectory  is  easily  plotted 
from  given  initial  conditions  on  the  phase-plane. 

b.  By  inspection  of  a  few  trajectories  the  sys- 

s tern's  behavior  is  determined.   For  example:  In  Fig.  2-2 

this  system  has  damped  oscillations  and  is  a  stable  system. 

Point  B  is  an  oscillation  max.  (0   =  1.2). 

o 

c.  In  this  example,  as  applied  to  the  case  of  a 
linear  servo  system  (C  is  constant),  isoclines  are  straight 
lines,  so  it  is  simple.   If  C  is  not  constant  more  labor 

is  required. 
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Fig.  2-2.   Phase-plane  plot  of  damped 

oscillation  obtained  by  iso- 
cline method  for  a  servomechanism. 
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B.    LIENARD  METHOD 

1.   Method  Description 

Lienard's  method  is  a  means  for  drawing  a  direc- 
tion field  in  the  phase-plane  to  guide  the  tracing  of  the 
trajectories  of  a  differential  equation. 

Using  the  isocline  method,  (see  Fig.  2-3) ,  draw  the 
zero-slope  isocline  on  the  phase-plane.   Let  P  be  any 
point  of  the  plane  at  which  the  slope  of  the  direction 
field  is  desired.   From  P  draw  the  line  parallel  to  the 
z-axis,  which  intersects  the  zero-slope  isocline  at  Q. 
Let  R  be  the  projection  of  Q  on  the  z-axis.   Now,  draw  a 
short  line  segment  through  P,  perpendicular  to  the  straight 
line  RP.   This  segment  is  tangent  to  the  direction  field 
of  the  integral  curves  of  certain  differential  equations. 

In  practice,  graphical  construction  of  the  entire 
direction  field  is  quite  easy  since  the  short  segments  at 
all  the  points  on  the  extension  of  the  line  PQ  can  be  drawn 
by  using  a  compass  centered  at  the  one  point  R.   Note  that 
this  construction  does,  indeed,  result  in  &  =   0  everywhere 
on  the  zero-slope  isocline. 

In  Fig.  2-3  the  two  angles  a.   and  (9  are  equal  because 
of  the  perpendicularity  of  their  sides  at  P  and  Q.   The 
tangent  of  each  angle  can  be  evaluated  in  terms  of  the  co- 
ordinates of  the  points  P,  Q  and  R  and  the  ratio  of  lengths 
of  unit  distance  along  the  two  axes.   This  length  ratio  is 
defined  by 
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Fig.  2-3.   Lienard's  construction  of 

phase-plane  direction  field. 
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_  unit  distance  for  v  ro_£\ 

unit  distance  for  z  * 


tan  oc  =  K  ~  =  KM  (2-7) 

dz 

where  M  is  the  slope  of  the  integral  curves  and/or  of  their 
direction  field 


tan^oR  =  il  ipt;  )  =  ~k^ —  (2_8) 


The  expression  Z_  =  f(VQ)  is  the  relation  defining  the 
zero-slope  isocline. 

Since  V  =  V  ,  P  is  any  point  of  the  phase-plane.  The 
subscript  can  be  dropped,  giving 


tan  P   =  *a=*  (2-9) 


since 

tan  &  =    tan  |9 


then 

tan  ct  =   KM  =  tan  (9  =  f£$  "Z 


KV 

M  =  f(^"Z  (2-10) 

K  V 

M  =  -=— 
dz 

2 
,,  dz    dv   dz    d  z  fo    -,  -■  v 

dt    dz   dt    -.,2 
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Combining  equation  (2-10)  and  (2-11) 

k2  ^t  - f  ( ft  ) +  z  -  °  <2-12> 

dt^       ac  " 

Any  servo  which  obeys  a  differential  equation  of  this 
form  has  a  phase-plane  representation  with  a  direction 
field  which  can  be  constructed  by  Lienard's  method. 

Note  that  the  coefficient  of  z  is  unity  and  that  K 
specifies  the  ratio  of  vertical-to-horizontal  unit  distan- 
ces on  the  phase-plane.   The  zero-slope  isocline  is  found 

d  z 

directly  from  the  differential  equation  by  setting  ~ 

dtZ 
equal  to  zero,  since  the  isocline  is  the  line  for  which 

d2z 
M  =  0.   So  equation  (2-11)  gives  o  =  0. 

at 

2.   Example 

Using  the  same  example  as  in  A,  a  linear  servo- 
mechanism  having  a  constant  damping  coefficient  C  =  0.25, 
and  natural  frequency  oj  =  10  radian/sec. 

d2fi         dG       « 

x9-   +  2cui     -rr^  +  ti>     8   =0 

,.2       n  dt     n   o 
dt 

Changing    to    the    required   form    (compared    to   Eq. (2-12)) 

.      d2e        0      de 

°    +    ^       TT^    +    0        =0 


2       ,L2  ui        dt  o 

U)  dt  n 

n 

d26  -,        d9 


100      dfc2  20      dt  o 


K   =  —  =   0.1 
n 
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,   dft   N      ~    d0       t    d0 
f  /^  o  \  _    2c_  o  _   _1_  o  _    V_ 

V  dt  J  u)        dt       20   dt       20 

from  which  the  equation  for  the  zero-slope  isocline  is 

obtained. 

ft   -  -  -¥- 
6o  "    20 

(See  Figs.  (2-4)  and  (2-5)). 

3 .   Discussion 

a.  Any  system  which  has  a  differential  equation 
like  Eq.  (2-12)  can  use  Lienard's  method  to  construct  a 
phase-plane  direction  field.   After  the  direction  field  is 
constructed,  one  can  draw  any  phase  trajectory  from  a  given 
initial  condition  on  the  phase-plane. 

b.  As  compared  with  the  isocline  method  in  Sec- 
tion A.   If  the  differential  equation  can  be  manipulated 
to  (2-^12)  form,  a  lot  of  labor  may  be  saved,  and  one  gets 
the  same  result.   It  seems  better  than  the  isocline  method 
in  this  particular  problem. 

c.  If  C  is  not  a  constant,  then  the  zero-slope 
isocline  will  not  be  a  straight  line. 

d.  Lienard's  method  (even  for  drawing  an  indivi- 
dual trajectory)  is  not  applicable  to  all  second-order 
non-linear  differential  equations. 
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V^0 


4    *. 


■Zero -slope 

isocline 
G0*-V/20 


Fig.  2-4.   Lienard's  construction 
of  the  phase-plane 
direction  field. 


\\\\\\\\V» 
\\\\\\\\N  N* 


Fig.    2-5.      Direction    field  obtained 
by   Lie'nard's   method. 


33 


C.    DELTA  METHOD 

1.   Method  Description 
Given 

x  +  G(x,x, t)  =  0  (2-13) 

where  (x,  x,  t)  may  be  given  analytically  or  graphically 
and  may  be  a  very  complicated  function. 

Rewrite  Eq. (2-13)  to  the  standard  delta  form.   This  is 

2 
done  by  adding  and  subtracting  the  term  p  x  from  Eq.  (2-13) 

2n     2 

x  +  [G(x,  x,  t)  -  p  x]  +  p  x  =  0 


Let 


2  2 

p  6  =  [G(x,  x,  t)  ,  p  x] 


so  that  the  linearity  difference  6  becomes 

6  =  \   [G(x,  x,  t)  -  p2x]  (2-14) 

P 


Eq.  (2-14)  will  make  (2-13)  appear  in  the  standard  fi  form 

x  +  p2(x  +  6)  =0  (2-15) 

The  geometric  significance  of  (2-15)  in  a  phase-plane 
representation  is  seen  by  introducing  the  phase-plane  co- 
ordinates 

x  =  x    and   —  =  v. 
P 


Then 


dv  dv    .  2      dv 

X   =  Pdt   =   PdIX   =   PVdx- 
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The  standard  second-order  equation  (2-15)  then  becomes 
a  first-order  differential  equation  in  x,  v,  and  fi,  namely 


v|+x+6=0  (2-16) 


^Z  =    _  Ix±6Jl     •  (2_17) 

dx  v 


Fig.  2-6  shows  in  this  case  that  the  slope  of  the  line 

x+6 
QP  will  be  ;   consequently  the  equation  of  the  normal 

to  QP  will  be  identical  with  Eq.  (2-17) .   The  instantaneous 
phase  trajectory  is  therefore  a  circular  arc  passing  through 
P  and  with  its  center  at  x  =  -6  or  at  Q. 

The  next  problem  is  therefore  to  relate  the  circular- 
arc  segment  beginning  at  P  to  the  finite  time  interval  At 


Since 


x 

v   =  — 


then 


dx 
dt 


P 


=  pv 


dt   =  —  (2-18) 

pv 

Fig.  2-7.   Let  the  average  radius  of  curvature  of  the  step 
be  given  by  QP  =  p. 
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Normol  to  OP 
P 


Fig.  2-6.   Exact  relationship, 


Fig.  2-7.   Stepwise  relationship 
for  constant  p. 


36 


Then 


d9  =  % 


-i 


2      2 

dx   +  dv 

v2  +(x+6)2 


dx 
v 


.*.  dt  =  -  d©  (2-19) 


Consequently,  an  integration  will  give 

t  -  t    =  ^(9-e  ) 

o   p    o 
pAt  =  A0  (2-20) 

Eq.  (2-20)  shows,  in  regard  to  the  angle  subtended  by  the 
circular  arc  described  by  the  phase  point  during  the  finite 
time  interval  At  . 

A  graphical  integration  of  the  differential  equation 
in  the  standard  delta  form  can  then  be  effected  by  a 
successive  construction  of  circular  arc  segments  whose 
centers  are  located  on  the  x-axis  at  various  step  value  of 
x  =  -6  (see  Fig.  2-8a) . 
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(a)  Change  with  positive  6.     (b)  Change  with  negative  fi. 


Fig.  2-8a. 
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If  x  had  been  taken  horizontal  and  positive  to  the 
right  with  v  positive  vertically,  positive  time  would  have 
been  represented  by  a  clockwise  angular  variation.   (see 
Fig.  2-8b) . 

The  definition  of  6  given  by  Eq.  (2-14)  transforms  the 
original  Eq.  (2-13)  into  the  standard  delta  form.   Curves 
representing  6  in  terms  of  the  phase  coordinates  x  and/or 
—  will  be  superimposed  on  the  phase-plane  so  as  to  permit 
choice  of  a  suitable  average  value  of  6  for  the  construction 
of  the  circular-arc  segments  making  up  the  phase  trajectory. 

2.   Example  (II-B-2) ,    Mass  on  Non-linear  Spring 

The  motion  of  a  constant  mass  attached  to  a  non- 
linear spring  which  becomes  relatively  stiffer  with  increas- 
ed deflection  is  described  by  an  equation  of  the  form 


^-~  +  25(1  +  0.1  x2)x  =  0 


Find  a  phase-plane  solution  curve  for  this  equation  by 
6  method,  with  the  initial  condition  that  at  t  =  0,  x  =  3 

and  -rr  =  0. 

dt 

Comparing  to  the  standard  delta-form  (Eq. (2-15) ) 


p2(x+6)  =  25(1+0.1  x2)x 


2     OC 

p   =25 
6  =  O.lx3 
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(a)  Change  with  positive  6.    (b)  Change  with  negative  6, 


Fig.  2 -8b. 


40 


First  draw  the  §  curve.   Phase  trajectory  construction 
proceeds  as  shown  in  the  figure.   The  initial  point  is 
x(0)  =  3  and  v(0)  =  0,  as  the  first  step.   It  is  assumed 
that  x  decreases  to  the  value  2.8  or  Ax   =  -0.2  as  shown. 
The  average  value  of  x  during  this  interval  is  x     =2.9, 
for  which  the  average  value  of  6  can  be  read  directly  as 
6     =2.4  approximately.   Thus  the  center  of  the  first 

circular  arc  is  located  at  point  x  =  -6     =  -24  and  v=  0; 

*  av       • 

the  radius  is  R  =  5.4  and  with  this  radius  an  arc  is  drawn 
to  the  point  where  x  =  2.8,  indicated  by  the  small  dot 
marked  x(0)  +  Ax     in  the  figure.   By  continuing  in  this 
manner,  the  entire  solution  curve  can  be  built  up  as  a  se- 
quence of  circular  arcs.  (Since  the  original  equation  in- 
volves no  damping,  the  solution  is  a  closed  curve  repre- 
senting a  periodic  oscillation.   (See  Fig.  2-9) . 
3.   Discussion 

a.  The  delta  method  is  limited  to  construction  of 
phase  trajectories  for  second-order  differential  equations. 

b.  In  starting  a  solution  by  the  isocline  method, 
the  entire  plane  must  be  fitted  with  line  segments  fixing 
the  slope  of  a  solution  curve.   If  only  a  single  solution 
curve  is  needed,  only  a  few  of  these  line  segments  are 
actually  put  to  use.   Considerable  simplification  would 
result  if  only  information  related  directly  to  the  desired 
solution  curve  were  used. 

The  delta  (6)  method  leads  more  directly  to  the 
desired  solution. 
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Fig.  2-9.   Phase-plane  diagram  constructed  by  fi  method, 
The  curve  for  non-linear  function  fi(x)  and 
the  final  solution  curve  are  shown.   Dots 
indicate  junction  of  successive  arcs  used 
in  building  up  solution. 
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c.   The  construction  of  phase  trajectories  will  be 
carried  out  for  various  magnitudes  of  steps.   So  the  accura- 
cy may  depend  on  different  amounts  of  constructional  work. 
4.   Notes 

a.   6  may  be  either  positive  or  negative.   The 
following  6  curves  for  five  systems  may  help  to  understand 
how  plots  of  curves  may  be  superposed  on  the  phase-plane. 

(1)  The  physical  pendulum   (See  Fig.  2-10) 

2 
x  +  p  sin  x  =  0 

■I  o 

6   =  — j  [G(x)-p  x]  =  sin  x  -  x 
P 

(2)  A   system  with  restoration  proportional  to 
the  cube  of  the  displacement   (See  Fig.  2-11) 

x  +  0  x3  =  0 

o  =  — p  X  -  X 

p 

(3)  A    linear    restoration   system  with  quadratic 
velocity   damping       (See    Fig.    2-12) 


x   +   q|x|x+px   =   0 

6   =   — 5-    [G(x,x)     -   p   x]=-^r   |x|x=q|^|   ^ 
P  P  P      P 
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2-10. 

>.  *"(sin;r-jr)  • 

Fig. 

The 

/ 

/I 

operative 

p 

displacement 

physical  pendulum  plotted  in  the  phase- 
plane.   Centers  0   and  0  correspond 
to  points  P  and  Q. 


.— -~    y 


Fig.  2-11.   The  operative  displacement  b   for  a 

system  with  restoration  proportional 
to  the  cube  of  the  displacement. 
Centers  0  and  C"  correspond  to 
points  P  and  Q. 
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(4)  A    system  with  positive    linear   and  with 
negative   quadratic    restoration   in  x      (See   Fig.    2-13) 

x   +  px-a|x|x  =   0 

6  =  — j(p  x   -   a    I  x   I  x   -  p  x)=   -  —~    I x   I x 
P  P 

(5)  A  system  acted  upon  by  nth -power  velocity 
damping  and  restored  by  first  and  mth-power  displacement  re- 
storation  (See  Fig.  2-14) 

••  ■  «1  -n-1  I  •     2     a  I  m-1  I     _ 

x  +  ff  |x     |x  +  px+  pri^x.  |  x  =  0 

2-n  vp     p     2  ' 

P     ^     *        P 


=  61  +  52 


(b)  The  b   curve  of  a  system  indicates  the  system's 
type  and  degree  of  non-linearity.   Thus,  if  a  system  is 
linear,  5  will  be  zero  unless  an  artificial  parameter  has 
been  introduced  deliberately. 

D.    SZEGO  METHOD 

1.   Method  Description 

Consider  the  second-order  system: 


(2-21) 
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Fig.  2-12.   The  operative  displacement  6  for 
a  system  with  quadratic  velocity 
damping.   Centers  0   and  0 
correspond  to  points  P  and  Q. 


I 


Fig.  2-13 


The  operative  displacement  6  for 
a  system  with  positive  linear  and 
negative  quadra toe  restoration  in 
displacement.   Center  0  corre- 
sponds to  point  P. 
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Fig.  2-14.  The  two  operative  displacements 
6,  and  fi~  for  a  system  with 
damping  proportional  to  the  nth 
power  of  the  velocity  and  with 
restoration  proportional  to  the 
first  and  the  mth  power  of  the 
displacement.   Center  0o  corre- 
sponds to  point  Q. 
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It  is  desired  to  identify  the  slope  of  the  trajectory 
of  this  system  at  every  point  of  the  phase  space.   We  con- 
sider the  topographical  system;   that  is  a  system  of  non- 
intersecting  differential  curves,  which  are  symmetric  with 
respect  to  the  origin. 

v  =  0(x,y)  (2-22) 

In  connection  with  the  function  (2-22),  construct  the 
curves  defined  by  the  equation 


&  -  I2  &  +  ¥     f£  -  °  (2-23) 

dt    Bx   dt    by   dt  x     ' 


and  Eq.  (2-21) . 

These  curves  we  called  contact  curves.   Substituting  Eq. 

(2-21)  into  (2-23) 

r>v       dy 

_|x_  =  _  _^t_   _   f(xfy)  (     ) 

dv        dx         y  v 

By         dt 

At  every  point  of  the  curve  (2-24) ,  the  corresponding  tra- 
jectory of  the  system  (2-21)  has  the  same  slope  as  the  curve 
of  the  family  (2-22)  at  the  same  point;    furthermore,  the 
curves,  representing  in  the  phase-plane  the  solutions  with 
odd-order  of  multiplicity  of  Eq.  (2-23)  divide  the  phase- 
plane  into  regions  in  which  -rr   has  fixed  sign.   On  the  con- 
trary, -rr     has  the  same  sign  on  both  sides  of  the  curves, 
representing  the  solution  with  even  order  of  multiplicity 
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of  Eq.  (2-23) .   Of  course,  in  the  latter  case  the  slope  of 
the  trajectories  at  every  point  on  these  curves  is  the  same 
as  the  slope  of  the  family  (2-22) . 

In  regions  where  -rr     is  negative,  the  representative 
point  of  the  system  (2-21)  describes  a  trajectory  in  the 
phase-plane  which  crosses  the  curves  of  the  family  curves 
(2-22)  in  the  direction  of  decreasing  v. 

In  the  regions  in  which  -tt  is  positive  the  trajec- 
tories cross  the  curves  of  the  family  (2-22)  in  the  direc- 
tion of  increasing  v. 

Choosing  different  families  of  curves  (2-22)  and  plot- 
ting the  corresponding  curves  -rr  =0   in  the  phase-plane, 
provides  enough  information  to  sketch  the  trajectories. 

If  the  scalar  function,  v  =  0(x,y)  changes  sign  in  the 
phase-plane,  i.e.,  if  the  equation  v  =  0(x,y)  has  some  real 
non-zero  solution  with  odd-order  of  multiplicity,  then  the 
curves  representing  these  solutions  must  also  be  considered 
as  boundaries  between  regions  with  different  sign  of  -rr  • 
This  must  be  done  in  order  to  take  into  account  the  varia- 
tion  in  sign  of  -jr   ,  due  to  variation  in  the  sign  of 
the  curves  of  the  topographical  system  (2-22) .   Notice  that 
the  trajectories  of  the  system  (2-21)  on  the  curves  v  =  0 
have  not  the  same  slope  as  curves  of  the  family  (2-22) . 

Consider  the  particular  topographical  systems: 

v  =  x  (2-25) 
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for 

which 

dv 

dt  =  Y 

and 

y1   =  y 

for 

which 

dv, 

dt   "  f(x'y) 

(2-26) 


(2-27) 


(2-28) 


The  subsequent  families  of  curves  are  now  considered  as 
topographical  systems. 


v2  =  x2  +  y2  (2-29) 

v3  =  x2  -  y2  (2-30) 

v4  =  xy  (2-31) 

Using  v,  =  y 
when   v,  =  0 

for     y  =  0  (2-32) 


and  from 

dvl   **l     dx   Svl  dy   dy 

dt     dx   dt    Sy   dt    dt 

then 

f(x,y)  =  0  (2-33) 
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2     2 

When  considering   v„  =  x   +  y  ;   it  is  seen  that   v_  =  0 

for   x  =  y  =  0   and 

dv 

=  0 


dt 

from 

dv 


2    n   dx  ,  _   dy 
"  =  2x  —  +  2y  - 


dt        dt     J    dt 

=  2xy  +  2yf(x,y)  =0 

then 

x  +  f(x,y)  =  0 

and 

y  =  0  (2-34) 

For 

2     2 
v3  =  x   -  y 

v   =  0    for   x  =  y 

and    x  =  -y  (2-35) 

Also 


dt 
But 


dv_ 

=  0  , 


dv„    Bv0   ,  av_   , 
3  ..  3   dx  _  3   dy 

dt     Bx    dt  By    dt 

„   dx  «   dy 

dt  J    dt 

=  2xy  -  2yf(x,y) 

=  0 
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Hence 

X 

= 

f  (x, 

-y) 

and 

Y 

= 

0 

(2-36) 


Finally,  considering  v.  =  xy,  gives 

v4  =  0 
for   x  =  0    and   y  =  0  (2-37) 

and 

dv 

=  0  , 


dt 

But 

dv.    dv         Sv.  , 
4  .  4   dx   4  dy_ 

dt   ~  dx    dt    By   dt 

dx  ,    dy 
J    dt      dt 

=  y  •  y  +  x  f (x,y) 

=  0 

Hence 

y2  =  -xf(x,y)  (2-38) 

2,      Example  (II-D-2) 

Consider  Van  der  Pol's  equation: 

x  =  y 

y  =  f(l-x2)y  -  x  (2-39) 
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dv. 

(a)  —  =   0     for    y  =   X   2 

1  -  x 

dv 

(b)  ^  =0     for    y  =  0 

x  =  ±  1 

dv2 
There  is  no  variation  in  the  sign  of  ,.    between  the  two 

2 

sides  of  the  curve  y   =0. 

dv., 

(c)  -r—   =0     for     y  =  0 


dt 


2x 
1  -  x 

(d)   ^=0     for    y  =  |  [(x2-l)iy(x2-l)2+  4] 

The  next  step  is  to  plot  the  curves  (a) ,  (b) , (c)  and 
(d),  (see  Figs.  2-15,  2-16,  2-17  and  2-18). 

In  order  to  use  the  method,  it  is  necessary  to  super- 
impose the  curves  drawn  in  Figs.  2-15,  2-16,  2-17,  and  2-18 

as  is  shown  in  Fig.  2-19. 

dv. 
The  best  procedure  now  is  to  check  the  sign  of   ,. 

(v  =  1,2,3,4)  numerically  at  one  point  of  the  phase-plane 

dv. 
and  utilize  the  result  for  identifying  the  sign  of   ,. 

(i  =  1,2,3,4)  in  the  different  regions  in  which  curves 

dv. 

-rr—  =   0   and  v  =  0(i  =  1,2,3,4)  divide  the  phase-plane. 

Since 

0  <  v      <   +    • 

_oo  <  v    <  +  °° 
1 

-00  <  V,  <  +  °° 

-«  <  v4  <  +  °° 

dv,    dv.. 

It  is  necessary,  before  checking  the  sign  of   -,,   ,  -.,   , 

dv4  ac     ac 

and  ,     at  one  point  of  the  phase-plane,  to  identify  the 

sign  of  v, ,  v_  and  v.    at  that  point. 
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v,«o 


H 


M 


v      =   0      for   y   =   0 


dv 


dt 


1         n  £  X 

— =  0      for  


1-x' 


generate    regions 
M,N,T,Q,R   and   S 


Fig.    2-15.      Curves      v,    =   0, 


and 


dv. 

dF 
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A, 


dt 
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H 1 1 1 1 h- 


+  1 


4 

dt 


Fig.    2-16 


dv. 

Curves    v»=    0    and  -rH 
2  dt 


=    0 


dv, 
dT 

y4 

x 


=   0      for 

=    0 
=    ±   1 


generate   regions 
a,  fi  and   v 
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B 
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:\ 


"3-0 


at 


3«o 


dv. 

Fig.    2-17.    Curves   v   =    0    and  -ttt 

J  dt 


=    0. 


V3 

dv. 
dt~ 


=    0      for      x   =   y 


=   0      for      y   = 


2x 
l-x: 


generate  regions 
A,B,C,D,E,F,G,H, 
I,j,  and  K 
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v,  =  0 


for  x 
dv 


0  ,  y  =  0 


dt 

for  y 


=  f  |_(x2-l)iv,(K2-l^  +  4 


generate  regions 
a,b>cd,e,f,g,  and  h 


=    0. 


VV$3a0 


*-X 


I»4=0 


J  -«•<?    \ 


m  ° 

£2=o    VOS?80! 


v3«o 


Fig.    2-19.       Example    (II-D-2) 
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Some  information  about  the  slope  of  the  trajectories 
at  every  point  of  the  phase-plane  is  thus  provided.   For 
example,  consider  the  point  P(-3;  1.5)  in  Fig.  2-20), 
dv. 


from 


from 


from 


dt 


dv. 
dF 


dv 

i 

dTF 


<  0  ,  requires  that  the  trajectory  move  into 

the  region  A'F'D'B* 

>  0  ,  the  trajectory  moves  from  P  into  the 
region  E'C'A'F' 

<  0  ,  requires  that  the  trajectory  move  into 


the  region  C'A'F'D'. 

These  three  conditions,  taken  simultaneously,  define 

the  angle  A'PF*  in  which  the  trajectory  which  originates 

at  P  must  lie. 

Since  the  slope  of  the  trajectories  is  known  exactly 

dv. 
at  the  points  on  the  curves   ;.   =  0  ,  it  is  possible  to 

draw  the  trajectories  starting  from  P  with  good  accuracy. 

(See  Fig.  2-19) 

Starting  for  example  from  point  #(-3,1.5)  contained 

dv. 
in  the  region  of  fixed  sign  of   ,    ,  namely  T,  a,    I  and  G, 

(See  Fig.  2-19) 

dv. 


dt 


<  0 


dv. 
dF 


<  0 
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dv. 
dF 


>   0 


dv 
dF 


<  0 


NOTE 


dv. 
dF 


<   0      in      T,    Q,    M   and   S 


dv. 
dF 


>  0      in      N  and  R 


dv. 
dF 


<  0      in      oi  and  y 


dv. 
dF 


>  0      in      j9 


dv. 
dF 


<   0      in      B,    D,    F,    H,    J   and  L 


dv. 
dF 


>  0      in      A,    C,    E,    G,    I    and  K 


dv, 
dt 


>  0      in     b,    d,    f   and  k 


dv, 

dt 


<   0      in      a,    c,    e   and  g 


3.   Discussion 

a.   This  method  is  a  simple  and  useful  method  for 
drawing  trajectories.   Also  it  is  essentially  a  generaliza- 
tion of  the  isocline  method. 
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b.  Accuracy  can  be  improved  by  considering  further 
topographical  systems. 

c.  It  is  worthwhile  to  notice  that  the  isocline 
method  is  a  particular  case  of  this  method.   And  it  corre- 
sponds to  the  particular  topographical  system: 


v  =  ax  +  a2y 


(-•;i».  J!+»;isll2) 


NOTE: 


dv  _  jVv   dx   _>v   dy_ 
dt  ;   ^x   dt    ^y   dt 


_  a   dx       dy_ 
"  al  dt  +  a2  dt 


=  0 

dy_ 
dy  ::   dt   .   _  fl 
dx    dx      a„ 

dt 


E.    MURTHY'S  NEW  APPROACH  TO  THE  PLOTTING  OF  PHASE-PLANE 
TRAJECTORIES 
1.   Method-Description 

This  method  is  a  new  and  simple  approach  for  the 
isocline  plotting  of  second-order  time-implicit  non-linear 
differential  equations  given  by  x  =  F(x,x'). 

a.   Isoclines  of  the  equation  x"+f (x ' ) +g (x)  =  0: 

x"  +  f(x')  +  g(x)  =  0  (2-40a) 
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Substituting  Mx '  =  x"  into  (2-24a) 


dx ' 
where     M  =  -= ;   gives 


Mx'  +  f(x')  +  g(x)  =0  (2-4Qb) 

(1)  Plot  from  the  given  differential  equation 
(2-40a)  the  curves  of  x"  vs  x '  on  the  x"x'  plane  for  con- 
stant (different)  values  of  x.   This  is  a  simple  matter 
since  it  amounts  to  shifting  the  curves  of  x"  =  f(x') 
vertically  for  different  values  of  x. 

(2)  Draw  a  number  of  straight  lines  x"  =  Mx' 
through  the  origin  (corresponding  to  the  different  values 
pf  M)  on  the  x"x'  plane;   this  can  be  done  quickly. 

(3)  With  any  particular  straight  line  in 
view,  say  for  M  =  M, ,  determine  all  possible  points  of 
intersections  of  this  straight  line  with  the  curves  obtain- 
ed from  (a).   Thus  a  number  of  values  of  x  and  x'  can  be 
obtained  corresponding  to  M  =  M, ,  which  when  plotted  on  the 
x'x  plane  and  joined  together  represent  the  isocline  for 

M  =  M,  .   The  procedure  is  repeated  for  M  =  M„,  M_.  ...  until 
a  fair  number  of  isocline  curves  are  obtained  on  the  phase 
plane.   The  respective  slope  lines  are  marked  on  these  iso- 
clines.  Joining  these  slope  lines  gives  the  different  tra- 
jectories or  the  entire  phase  portrait  for  Eq.  (2-40a) . 
(b)  Isoclines  for  the  equation 

x"  +  f(x')h(x)  +  g(x)  =  0  (2-41) 
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The  above  procedure  can  still  be  applied  except  that  the 
curves  of  x"  =  f(x')  are  now  multiplied  by  a  scale  factor 
[due  to  the  h(x)  term]  corresponding  to  different  values 
of  x. 

(c)  Notes 

(1)  This  method  is  more  convenient  in  deter- 
mining the  solution  curves  for  the  following  differential 
equations : 

x"+  x'h(x)  +  g(x)  =  0 

x"  +  x'     +  g(x)  =  0  (2-42) 

x"  +  xx'    +  g(x)  =  0 

when  compared  with  existing  methods. 

(2)  This  method  can  also  be  used  for  the 
following  nonlinear  differential  equations: 

h(x)x"  +  K(x)  f(x')  +  g(x)  =  0 

h(x')x"+  K(x)  f(x')  +  g(x)  =  0  (2-43) 

x"  =  F(x,  x') 

(3)  The  Lienard  method  is  a  step-by-step  con- 
struction for  drawing  any  particular  trajectory  correspond- 
ing to  a  given  set  of  initial  conditions.   The  present 
method  determines  the  entire  phase  portrait  indicating  the 
behavior  for  all  initial  conditions. 

(4)  Variations  of  the  nonlinear  restoring- 
force  term  g(x)  can  easily  be  studied. 
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(5)  The  method  is  independent  of  the  nature 
of  the  functions  f,g,h,k,  i.e.,  whether  logarithmic,  alge- 
braic, or  transcendental. 

(6)  For  nonlinear  differential  equations  of 
the  type  x"  =  f(x),  the  x"-x'  curves  become  horizontal 
straight  lines  for  constant  values  of  x.   The  intersections 
of  these  straight  lines  with  x"  =  Mx '  lines  give  the  points 
on  the  phase  plane  for  any  chosen  isocline  slope. 

(7)  The  isoclines  corresponding  to  M  =  0  and 
M  =  °°  can  generally  be  plotted  quickly. 

(8)  Sometimes  an  x"-x'  curve  (for  a  constant 
x)  and  an  x"  =  Mx'  line  may  run  so  close  to  one  another 
that  the  isocline  curve  on  the  phase  plane  is  nearly  ver- 
tical in  that  range  of  x'. 

2.   Example 

Consider  a  system  with  nonlinear  damping  and  re- 

2     2 

storing  forces  described  by  x"  +  x1   +  x  =0. 

(a)  Draw  x"-x'  curves  (see  Fig.  2-22). 

(b)  Draw  Mx '  =  x'  corresponding  to  various 
values  of  M  (a  number  of  straight  lines  through  the  origin. 

(See  Fig.  2-22) 

(c)  Draw  isoclines  and  phase  trajectories  for 

2     2 
x"  +  x'   +  x   =  0. 
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Fiq.  2-22.   x"A'  curves  and  x"  =  Mx '  lines  for 

2    2 
the  equation  x"  +  x*  +  x   =0. 
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Considering  the  line  for,  say  M  =  +  5,  the  points  of  in- 
tersections are  determined  as 


(x   =   0      , 

x' 

=      0) 

(x   =    2       , 

X  ' 

=    -1) 

(x  =  JT  , 

x' 

=    -7) 

(x  =  Js   t 

X' 

=    -1.4) 

(x   =  Js    , 

X  ' 

=    -2.0) 

(x   =  J6    , 

x' 

=    -2.5) 

These  points  are  marked  on  the  phase-plane,  joined  and  the 
slope  lines  are  marked  on  it  (Fig.  2-23) .   Similarly,  dif- 
ferent isoclines  are  obtained  for  different  values  of  M. 
The  trajectories  are  then  drawn  as  shown  in  Fig.  2-23. 
3.   Discussion 

(a)   This  is  a  simple,  all  graphical,  useful  and 
general  method  to  obtain  the  phase-plane  trajectories  cor- 
responding to  all  pertinent  initial  conditions  of  second- 
order  autonomous  systems  (which  are  not  amenable  to  the 
usual  isocline  method)  by  making  use  of  the  isocline 
method.   It  reduces  the  labor  of  calculating  isoclines  for 
many  non-linear  differential  equations  and  thus  makes  the 
isocline  method  a  much  more  useful  tool. 
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V  ■+*+++ >+kH^co 


Fig.  2-23.   Isoclines  and  phase  trajec- 

2     2 

tones  for  x"  +  x'   +  x  =  0, 
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F.    THE  ACCELERATION-PLANE  METHOD 

1.  Me  thod-De  sc  r ip  tion 

It  is  based  on  the  fundamental  equation 

^  =  ~   ,  (2-44) 

dx   v 

where  w  =  w(t,  v,  x)  or  the  acceleration  is  a  function  of 
t,  v  and  x.   For  the  sake  of  simplicity,  assume  that  a  con- 
stant driving  force  is  applied. 

Thus      dv  =  w(v,x)  (2-45) 

dx     v 

The  method  calls  for  preliminary  plots  of  w  vs  x  for 
constant  values  of  v.   This  auxiliary  w-x  plot  for  constant 
values  of  v  is  not  a  trajectory.   It  characterizes  the  sys- 
tem just  as  the  original  differential  equation  character- 
izes the  physical  condition  of  the  system. 

This  method  is  applicable  to  both  linear  and  non-linear 
problems.   In  the  w-x  plot,  though  constant  driving  func- 
tions can  be  included,  the  time-varying  forcing  or  driving 
functions  must  be  added  for  different  values  of  time. 

2.  Example    x"  +  0 . 5x '  +  x  =  0. 

For  v  =  x'  =0  ,  then  x"  =  -x  ,  and  the  x"  -  x 

plot  or  w-x  plot  is  a  straight  line  with  a  slope  -1  . 

For  v  =  x*  =1,  then  x"  =  -x  -0.5,  so  this  straight  line  is 
displaced  from  the  first  line  by  -0.5. 

In  Fig.  2-24 (a)  six  points  of  construction  are  shown 

for  the  v-x  trajectory,  starting  from  P  where  t  =  0, 

x   =2.5  and  v   =  0 . 
o  o 
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Taking  x,  as  the  center,  and  starting  from  w. ,  an  arc 
can  be  drawn  in  the  clockwise  direction  so  that  the  arc 
strikes  R,  on  the  horizontal  axis.  From  R,  draw  a  line  to 
P,  where  x  =  x  and  v  =  v. ;  the  perpendicular  line  to 
P-.R-,  gives  the  correct  slope  for  the  v-x  trajectory  at  P..  . 
If  the  perpendicular  line  gives  a  slope  that  agrees  fairly 
well  with  the  assumed  slope,  point  P,  can  be  considered  as 

the  desired  point  which  falls  on  the  v-x  trajectory. 

Wl 
Slope  of  perpendicular  line  = 


Vl 
If  the  assumed  slope  does  not  agree  with  the  slope  of 

the  perpendicular  line  obtained  from  the  above  construction, 
a  new  slope  should  be  assumed  and  the  above  procedure  re- 
peated until  the  assumed  slope  agrees  with  the  slope  of  the 
perpendicular  line  (see  Fig.  2-25).   On  Fig.  2-24 (a),  P. 
has  the  value  x  =  x,  and  v  =  v,  and  v,  is  chosen  as  -0.5. 
Thus  if  one  prolongs  the  vertical  line  from  x..  downward 
through  P,,  this  line  should  intersect  with  the  v  =  -(h) 
line  in  w-x  plot  (shown  in  dot  lines)  at  a  point  w..  .   By 
connecting  the  points  w  ,  w.  ,  w_  ...  etc.,  the  w-x  trajec- 
tory is  constructed  at  the  same  time  (see  Fig.  2-24(b)). 
3.   Discussion 

This  method  is  applicable  to  both  linear  and  non- 
linear problems.   One  draws  the  w-x  trajectory  first,  then 
gets  the  v-x  phase-plane  trajectory.   Each  time  only  one 
trajectory  is  obtained,  so  this  is  a  useful  method. 
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Fig.    2-25.      Acceleration-plane  method, 
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G.    DEEKSHATULU ' S  METHOD 

This  method  suggests  the  use  of  simple  transformations 
for  second-order  non-linear  differential  equations  to  effect 
rapid  plotting  of  the  phase-plane  trajectories. 
1.   Method-Description 

a.   Transformations  for  second-order  autonomous 
systems 

A  general  second-order  non-linear  autonomous 
system  is  governed  by  a  differential  equation  of  the  form 

x  =  f(x,x)  (2-46) 


Let 


then 


dx 
y  =   dt 


2 

d  x  _  dy_  _  dy_   dx 

,,2    dt    dx    dt 
dt 


=   M 


dx 
dt 


Change  Eq.  (2-46)  to 

Mx  =  f(x,x)  (2-47) 

To  draw  the  isocline  for  an  assumed  value  of  M,  equal 
to  M, ,  say,  x  has  to  be  found  as  a  function  of  x  and  M,  and 
this  in  many  cases  is  not  possible;   thus  the  only  attempt 
made  is  to  find  the  points  of  intersections  of  such  an 
isocline  with  a  given  set  of  straight  lines  or  with  simple 
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standard  curves.   Consequently,  on  the  phase-plane,  if  one 
considers  the  intersections  of  the  isocline  for  slope  M  =  M, 
and  the  set  of  straight  lines  x  =  Kx (K  =  ±  1,  K  =  ±  2,  ...), 
Eq.  (2-47)  becomes 

JV^Kx  =  f(x,Kx) 

or 

M2  =  F(x,K)  (2-48) 

For  a  specified  K,  the  solution  of  Eq.  (2-48)  (algebraic 
equation  with  constant  coefficients)  for  real  roots  of  x 
gives  the  points  of  intersection.   But  it  is  easier  to  find 
the  values  of  the  slope  M  (for  various  values  of  x)  on  the 
lines  x  =  Kx  than  it  is  to  try  to  obtain  the  isocline  curve 
for  a  given  value  of  the  slope.   Even  under  this  transforma- 
tion x  =  Kx,  every  quadrant  of  the  phase-plane  is  uniquely 
defined.   It  is  advisable  to  find  any  symmetry  present  in 
the  phase  portrait  for  a  given  differential  equation  since, 
if  it  exists,  it  reduces  labor  and  is  tested  as  follows: 

(1)  Symmetry  about  the  x-axis 

If  for  K  =  positive  (negative)  K,  and 
x  =  positive  (negative) x. ,  the  slope  M  =  M, ,  then  for 
K  =  negative  (positive)  K..  ,  and  x  =  positive  (negative)  x^, 
the  slope  M  =  -M, . 

(2)  Symmetry  about  the  x-axis 

If  for  K  =  positive  (negative) K,  and  x  = 
positive  (negative) x, ,  the  slope  M  =  M2 ,  then  for  K  =  nega- 
tive (positive)K, ,  and  x  =  negative  (positive) x_,  the  slope 
M  =  -M2. 
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The  procedure  for  determining  M,  on  the  phase-plane 
for  any  chosen  value  K.  and  x..  can  be  viewed  in  the  phase 
space  (with  x,  x,  and  x  as  coordinate  axes) .   Determine  the 
curve  of  intersection  of  the  surface  defined  by  Eq.  (2-46) 
with  the  plane  x  =  K,  x,  then  determine  the  simple  point  of 
intersection  P  of  this  curve  with  the  plane  x  =  x, .   The 
slope  of  the  plane  passing  through  the  x-axis  and  the  point 
P  with  x-x  plane  determines  the  value  of  M, .   Similarly, 
other  values  of  M  are  obtained  by  taking  different  values 
of  x  and  the  procedure  is  repeated  for  each  value  of  K. 

Other  transformations  such  as  x  =  Kx  ,  x  =  Kx  ,  x  or 
x  =  K,  etc.,  can  also  be  considered. 

This  procedure  is  general,  systematic,  and  comparative- 
ly simple  when  compared  with  the  method  of  determining  slope 
lines  arbitrarily  at  points  for  the  phase  portrait  construc- 
tion. 

b.   New  Planes  for  Certain  Second-order  Time- 
Varying  Systems. 

Time-varying  elements  occurring  in  physical 
systems;   for  example,  the  variable  mass  of  a  guided  missile 
can  be  though  of  as  function  of  time  since  its  weight  de- 
creases as  the  fuel  is  consumed. 

On  the  phase-plane,  solution  of  this  kind  of  problem 
makes  use  of  step-by-step  procedures  by  considering  simul- 
taneously the  x-x  and  x-t  (solution)  planes.   Professor 
Deekshatulu  suggests  some  new  planes  whose  coordinate  axes 
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are  explicit  functions  of  time  and  which  effect  an  easy 
and  rapid  solution  of  at  least  certain  classes  of  second- 
order  non-linear  time-varying  systems  not  easily  amenable 
otherwise  by  the  phase-plane  methods. 

c.   Suggestions  for  Solving  Certain  Higher-order 
Systems 

Assuming  non-linear  transformations  of  the 
form 

£  =  e  +  af(e)   or   t?  =  e  +  bf(e)+  ae    (2-49) 

third  and  fourth-order  time-varying  systems  described  by 
the  following  equations: 

t2e  +  \  +  f  (4)  =  o 

2~€i  +  \  +  f  U)  =  o 


(2-30) 


t  t?  +  tf?+  f  (n)  =  o 


\     2trj  +  f?  +  f(r?)    =  o 

can  first  be  solved  by  the  method  described  and  then  be 
solved  as  a  first-order  or  second-order  nonautonomous  sys- 
tem.  Note  that  the  equation  of  the  form  e  +  af(e)  =  h(t) 
can  be  solved  on  the  e-t  (solution)  plane. 
2.   Examples 

(Example  1)   Consider  a  non-linear  system  with  dis- 
sipation  proportional  to  x  ,  given  by  x  +  hx   +  x  =  0 
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Let 

dx 


then 


dt   =   Y 


2 

d  x  _   dy_  _  dy  dx  _    dx 

,,2     dt  :  dx   dt      dt 
dt 


Thus 


,  .2 
x  +  hx   +  x 


=  Mx  +  hx   +  x  =  0 


or     hx   +  Mx  +  x  =  0 


Both  roots  of  this  equation  are  valid.   On  the  other  hand, 
if  one  substitutes  x  =  Kx  this  equation  becomes 

h(Kx)2  +  MKx  +  x  =  0 
or 

-1-K2hx 

M  =  

K 

which  is  simpler  to  plot  for  any  given  values  of  h  and  K. 
Since  the  phase  trajectories  are  found  to  be  symmetrical 
about  the  x-axis,  only  the  upper  half  of  the  portrait  is 
shown  in  Fig.  2-26  for  the  case  h  =  1. 

(  «-^ 

for   h  =  1     < 

(   x  =  Kx 
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etc. 


(Example  2)   Let  it  be  required  to  find  the  behavior 
of  the  following  time-varying  system: 

t2x  +  tx  +  x3  =  0 

for  all  pertinent  initial  conditions  (with  t  4   0) . 

The  task  would  be  enormous  if  it  were  to  be  attacked 
by  the  phase-plane  method.  Using  the  tx-x  plane  the  dif- 
ferential equation  reduces  to 

3 
x 

tx  =  y  =  -  or 

Note: 
Let 

tx  =  y 

then 

d£   =   ^^   =   Mx  =  x  +  tx 
dt      dx   dt 

2       3 

Mtx  =  tx  +  t  x  =  -x 

or 

3 

tx  =  y  =  -  — 
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Taking  different  values  for  M  the  isoclines  (enlarged 
or  reduced  scale  plot  of  the  cubic  curve)  are  drawn  and  the 
respective  slopes  marked  on  them  with  their  direction  (see 
Fig.  2-27)  . 

3.   Discussion 



2 

a.  Simple  substitutions  like  x  =  Kx,  Kx  , 

o      o      o  f*  ■ 

x   +  x  —  y      in  the  isocline  equations  6f;  second-order  sys- 
tems simplify  plotting  of  the  phase  trajectories  to  a  great 
extent.   The  method  is  a  useful  supplement  of  the  isocline 
method.  ,  ,  .   .,. 

b.  Since  there  are  not  many  methods  to  evaluate 
quickly  the  bahavior  of  second-order  non-linear  time-vary- 
ing systems  for  various  pertinent  initial  conditions  and/or 
for  different-.- non+OLinea.'r  characteristics,  Deekshatulu 

-   suggests;  in- general  time-varying  coordinate  planes,  the 

.  2 
use  -of  tx-x  and  tx  ^x  planes, .f.p.35;j  obtaining  .easily  the  solu- 
tion of  certain  classes  of,  nonlinear*  time-varying  systems. 

r 

C 

H.    PELL'S  METHOD  -f-    ~ 

1.   Method-Description 

Quite  generally,  second-order  nonlinear  autonomous 
equations  of  the  form 


n 


d  x  ,  dxN  n 
—  2  +  9(*<  dt>  =  0 
dt 

.  f 

may  be  written  as 

'  +  <*  (■#)  +  f  (x)  =  0.  (2-31) 


,.2    w  vdt 

dt 


x     ...  , 

m 
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Let 


dx 
y  =  dt 


then 


Thus 


2 

d  x  _  dy_   dy_   dx 

,.2    dt  :  dx    dt 
dt 


2 
^-f  +  0    (~)  +  f  (x) 

dt        c 


dy    dx  ,  ./dx,  ,  _  ,  .    _ 

=   j   *  T7T  +  0("jr)  +  f(x)  =  0 

dx    dt      dt 


or 


-<*($:)  -  f(x) 

dy  _    dt v  ' 

dx        dx 
dt 


til    ~    fM  (2-32) 

y 


Now,  as  shown  in  Fig.  2-28,  plot  -6(y)    versus  y  and 
-f(x)  versus  x.   Given  the  initial  conditions  x(0)  and  y(0), 
the  construction  of  a  segment  of  the  trajectory  through  the 
initial  condition  proceeds  as  follows: 

a.  Move  vertically  from  IC  and  find   f[x(0)]  ; 
lay  this  off  on  the  x-axis  from  x(0)  by  use  of  a  compass  or 
a  45°  triangle.   This  establishes  point  A. 

b.  Draw  line  AB,  where  OB  =  y(0) . 

c.  Find  -#[x(0)]  by  moving  horizontally  and  draw 
CD  parallel  to  AB. 
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d.   Now  the  hypotenuse  (D)  (IC)  of  a  right  triangle 
whose  base  is  the  distance  DE  =  #[y(0)]  +  f[x(0)]  and  whose 
altitude  is  the  distance  (E)  (IC)  =  y(0)  is  established. 
The  slope  of  the  hypotenuse  is 


(E)  (IC)   _  y(o)  f2_33> 

DE      "   0[y(O)3  +  f[x(0)]  U  JJ/ 


e.   The  negative  reciprocal  of  Eq.  (2-33)  is  the 
slope  required  in  Eq.  (2-22) ;   thus  a  line  perpendicular  to 
the  hypotenuse,  passing  through  the  initial  condition  forms 
a  segment  of  the  required  trajectory.   The  usual  precautions 
must  be  observed  as  to  the  length  of  the  extrapolated  seg- 
ment. 

Note:   The  major  advantages  of  this  method  are: 

(1)  Its  simplicity.   Only  two  triangles  are 
needed  to  carry  out  the  solution  once  the  two  functions  are 
plotted. 

(2)  It  eliminates  the  trial  and  error  of  the 
delta  method. 

2.   Example 

Consider  the  equation  (the  hydraulic  speed  control 
system) 

2 
^-f  +  Cl  |f  +  U   2(l+a2e2)e  +  S|  =  0         (A) 
dt2     1  dt     °  C2 


where 


c,  =  viscous  damping  coefficient 
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c„  =  inertial  coefficient  related  to  volume  of  nozzle 

chambers 
K  =  gain  of  amplifier 
Ui     =   natural  frequency  of  vibration  of  flapper  for 

small  amplitudes. 
The  equivalent  flapper  spring  is  assumed  to  be  linear 
Eq. (A)  may  be  rewritten  as 

2 

d  e  .     de  ,  ,,  2M  ,,2  2,     _  ,_. 

— 2    cl  dt     1  ^     e  )e  =  0  (B) 

dt 

Let   T  =  o>,  t   to  change  time  scale;   thus  Eq.  (B)  becomes 

2     c 
d  e  ,     de  ,  ,,  ,,2  2.  _  ,_. 

— 2    +  1Z   d?  +(1+b  e  )e  =  °  (C) 

dT      1 

Choose  the  following  more  or  less  arbitrary  set  of 
coefficients  designed  to  emphasize  the  construction  tech- 
nique : 

c,  =  164  U)     =   377 

1  o 

c2  =  0.01  a2  =  0.858 

K   =  1000 
on  a  consistent  set  of  units.   Then 


2 

%2  + 

K 
C2 

b2      = 

2 
o 

*i2 

0.58 

2 
a        = 

0.5 

fi 

0.33 
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Let 

x  =  e 


de 


Eq. (C)  becomes 


dy  ..  -0.33y  -(1+  0.58x2)x 
dx  v 


0(y)  and  f (x)  may  be  identified  and  plotted  and  the 
construction  completed  as  shown  in  Fig.  2-29. 
3.   Discussion 

This  method  is  satisfactory: 

a.  If  only  a  single  solution  curve  is  needed. 

b.  Only  a  few  of  these  line  segments  are  actually 
put  to  use. 

I.    SOMAYAJULU'S  METHOD 

In  this  method  a  graphical  procedure  is  suggested  for 
drawing  phase-plane  trajectories  making  use  of  tangent 
slopes  instead  of  normals  (i.e.,  the  negative  reciprocal  of 
tangent  slope  values)  as  is  done  now  in  the  Lienard  con- 
struction. 

1.      Method-Description 

A  straightforward  construction  is  presented  for 
more  general  cases  to  which  the  isocline  method  is  not 
applicable  and  construction  of  individual  trajectories  is 
required. 
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\/C[x  (0),  2.(0)] 


Fig.    2-28.       Pell's   Method   of  phase  .-trajectory 
construction. 


CB=AD 


Fig.  2-29. 


Example  construction  using  Pell's 
method  for  the  hard-spring  system. 
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a.   Let  the  slope  of  a  second-order  system  be 

given  by 

dx    x 

M  =  -z—   =  — 
dx    k 


f1(x)  +  f2(x) 


gx(x)  +  g2(x)  (2"34) 


Assume   f 1 ,  f„,  g,  ,  and  g„  are  all  continuous,  single-valued 
and  of  the  odd  characteristic  type;   these  curves  are  first 
drawn  on  the  phase-plane. 

Referring  to  Fig.  2-30,  the  initial  point  P(x  ,x  )  is 
marked.   Make  PQ  =  f ..  (x  ),  QR  =  f  ?  (x  ).   Similarly  the 
length  RST  is  made  equal  to  g, (x  )  +  g„ (x  ).   The  angle  PRT 
is  a  right  angle.   Note  that  the  horizontal  line  RST  is 
drawn  to  the  left  of  PR  since  the  slope  is  known  to  be  po- 
sitive from  the  nature  of  the  functions  f -.  ,  f~,  g.  ,  and  g~ 

and  the  initial  conditions  x   and  x  .   A  small  length  of  the 

o      o  J 

straight  line  joining  P  and  T  gives  the  required  slope  at  P. 
The  field  direction  is  marked  on  the  basis  that  with  x  posi- 
tive x  should  increase.   It  is  seen  in  the  procedure  that 
it  is  necessary  to  first  determine  the  direction  of  the 
normal  to  the  phase  trajectory. 
b.   Consider  equation 

x  +  f(x)  h(x)  +  g(x)  =  0  (2-35) 

This  equation  is  considered  to  be  of  practical  impor- 
tance because  of  the  presence  of  nonlinear  damping  and 
nonlinear  restoring  forces. 
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Let 


dx      , 

y  =  he    and 


rewrite  the  equation  (2-35)  in  form  (2-36) 


=  _  f(x)  h(x)  +  g(x)  (2-36) 


M 

x 


The  plots  of  f(x),  h(x)  and  g(x)  are  shown  in  Fig.  2-31. 

•  2        2 

Also  a  parabola   x    or   x   (only  one  half  will  suffice)  is 

drawn  in  addition  to  the  above  curves. 

Let  the  initial  point  be  at  P (x  ,  x  )  ; 

r  o   o 

then: 

OS  =  x      ,       OQ  =  x 
o  o 

QR  =  h(x  )  ,   SB  =  f  (x  ) 

**       o  o 

and   SA  =  g (x  ) 

3   o 

Make  SM  =  QR. 

Draw  a  circle  with  BM  as  a  diameter  and  let  it  intersect 

the  x-axis  at  the  points  C,  C. 

Then 


SC  =  Jf  (x  )  h(x  ) 
o     o 


With  the  help  of  the  parabola  the  squared  length  of  SC 
is  obtained  and  is  marked  as  PN. 

Make     NL  =  g(x  ) 
3  v  o 

and     LK  =  x 
o 
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The  segment  line  through  P  joining  P  and  K  represents  the 
required  tangent  slope  of  the  phase  trajectory  at  P.   The 
direction  of  motion  is  decided  to  be  downwards.   The  new 
values  x,  and  x,  at  the  end  of  this  segment  are  determined 
and  the  procedure  repeated. 

Note:   (a)   In  both  the  above  mentioned  methods  it  is 
necessary  to  take  equal  scales  on  the  x  and  x  axes. 

(b)   The  accuracy  in  plotting  is  enhanced  by 
considering  small  segments. 

2.   Example 

x  +  x   +  x  =  0 


with  initial  condition   x   =  -4   and  x   =  0 

o  o 


2 

Rewrite  x  +  x   +x=0   to  the  form 


-x  -  x    x+x 

M  =  =  

X        -x 


and  draw 

f1   =  x        f2  =  x 

q1    =  -x 

on  the  phase-plane  (see  Fig.  2-32)  (using  Somayajulu's 
method  as  in  I-l-a) . 
3 .   Discussion 

This  method  is  based  on  geometrical  relationships, 
so  a  precise  step-by-step  procedure  must  be  followed. 
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J.    THE  E-FUNCTION  METHOD 
1.   Method-Description 


^-|  +  f  (|*  ,  x)  +  g(x)  =  h(t) 

dt 


(2-45) 


Let  the  Kinetic  energy  and  potential  energy  be  defined  as 
[%)x2    and  Q(x) 

where 


«x 
Q(x)  =    g(x)dx. 
o 


(2-46) 


Define  the  total  energy  as 


E  =  -|(x)2  +  Q(x). 


(2-47) 


Differentiating  E  with  respect  to  time 


DE  _  dxo  d  x    dQ(x) 

dt    dt"  ,,2  dt 

dt 


(2-48) 


E  =  x  •  x  +  Q(x) 


Differentiating  E  with  respect  to  x 


dE   _.     d 
dx        dx 

L    2(X)      _ 

+   dQ(x) 
dx 

dx    [     d  fdx,   I       dQ(x) 
dt    L  dx^dt;   J         dx 

=   x   + 

g(x) 

(2-49) 
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From  Eq.  (2-49) 

2 

d  x    dE     .  . 

— 2  =  dx"  "  g(x) 

dt 

Substituting  this  into  Eq.  (2-45)  yields 


||  -  g(x)  +  f(|f  ,  x)  +  g(x)  =  h(t). 


or 


||  =  h(t)  -  f (x,  x)  (2-50) 


From  Eq.  (2-47) ,    an  E-x  plot  can  be  made  for  different 

values  of  x.   These  curves  may  be  called  equi-velocity 

curves,  as  x  =  v,  and  they  represent  different  kinetic 

energy  levels.   For  x  =  0,  the  E-x  curve  is  simply  E  =  Q(x) 

where  Q(x)  represents  the  potential  energy.   At  x  =  x, ,  the 

difference  between  any  two  E-x  curves,  say,  between  the 

curves  for  v  =  v,  and  v  =  v   is  equal  to  the  difference 

1  o     ^ 

(v,2-  vo2)/2. 

The  value  of  v  =  x   is,  from  Eq.  (2-47) 

x  =  v  =  ±  <v/2[E-Q(x)  ].  (2-51) 

Substituting  Eq.  (2-51)  into  (2-50), 

||  =  h(t)  -f["x,  +V2[E-Q(x)  ]]   for   v  >  0 


(2-52) 


||  =  h(t)  -f[x,  -72[E-Q(x)]]   for   v  <  0 
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The  E-f unction  method  is  thus  a  method  for  construct- 
ing the  E-x  trajectory  according  to  Eq.  (2-50) .   See  Fig. 
2-33,  the  slope  of  AB  is  equal  to  h(t)-f(x,v   ),  where 

Vl+  V2 

v    =  x 

av      2 

Xl+  X2  .   . 

and  x  should  be  corrected  to  ~ ,  even  though  it  is 

initially  taken  as  equal  to  x, . 

2.   Example:   Van  der  Pol  Equation 

o 
x-(l-x)x+x=0 

E-x  plot  (see  Fig.  2-34) . 

From  E-x  plot  get  the  v(v  =  x) 

Note:   As  compared  with  the  isoclines  and  trajectories 
of  the  Van  der  Pol  equation  obtained  from  the  digital  com- 
puter, errors  exist  in  this  method.  (See  Fig.  2-35) 

K.    THE  SLOPE-LINE  METHOD  (See  Fig.  2-36) 
Given  an  equation 

f  =  '« 

or   dy  =  f(x)dx,  the  incremental  relations  may  be  written 


Ay  =  f   Ax  (2-53) 

■*    av 

fl  +  f2 
where  f   = represents  the  average  value  of  f  (x) 

av      2       ^ 

during  the  interval. 

Thus 

A     Ax  -   ,  Ax  - 
Ay  =  —   fl  +  T   f2 

(2-54) 
=  AYl  +  Ay2 
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AO»l    CO«h(t)-|  (X  .-V5-) 


Fig.  2-33.   The  E-function  method, 


•**      -t       -t$     •!       S      a       -s       i        »s   «  a 


Fig.  2-34.   Example  (II-J-2) . 
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Ax 
Let  an  angle  6  be  chosen  such  that  tan  9   =  -~-  ,  for 

a  chosen  value  of  Ax,  then 

Ay, 

tan  0  =  —-±- 
rl 

=  -^  (2-55) 

r2 


Starting  from  an  ordinate  f. ,  make  an  angle  0  and  draw 
a  dotted  line  be  to  intersect  the  horizontal  line  at  c. 
Then   ac  =  Ay.   From  c,  draw  another  dotted  line,  making 
the  same  angle  0   with  the  vertical  line  erected  from  c, 
until  another  ordinate  f_  is  intersected.   As  shown  de  =  fp 
can  be  projected  from  the  f (x)  curve  on  the  left  of  the 
figure. 

Then  ce  =  Ay„,  and  ae  =  Ay,  corresponding  to  the  in- 
crement Ax  shown  on  the  left  of  the  figure. 

This  method  can  be  applied  to  the  solution  of  simulta- 
neous first-order  differential  equations,  linear  or  non- 
linear.  For  a  first-order  nonlinear  equation 


+  g(x)  =  h(x)  (2-56) 


dy 
dx 

the  incremental  relations  give 


Ay  =  [h(x)  -  g(y)   ]Ax  (2-57) 

d  V 
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isoclines 

and 

trajectories 


;     -I    Mill 

■*  -1 


Fig.    2-35. (Example    II-J-2) 


Fig.    2-36.      The    slope-line  method, 
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For  example,  in  a  problem  of  surge  tank  transients,  the 
equations  are  of  the  following  form: 

-  f|=  y ■  +  f (u)  (a) 


|*  =  v(T,  y)  -  u  (B) 


When  the  equations  are  in  normalized  form  both  the 
horizontal  and  vertical  systems  of  slope  lines  can  be 
drawn  with  the  same  slope.   Eq„ (A)  represents  acceleration. 
Eq. (B)  is  obtained  from  the  condition  of  continuity. 

L.    COMPARISON 

It  is  readily  seen  that  some  methods  are  easier  than 
others  but  may  be  of  limited  usefulness,  applying  only  to 
some  particular  problems.   Some  of  them  are  useful  but  com- 
plicated, so  that  it  is  easy  to  make  a  mistake.   It  seems 
the  most  useful  and  general  method  is  the  isocline  method, 
when  one  solves  problems  by  hand  labor.   Deekshatulu ' s  and 
Murthy's  methods  are  useful  supplements.   The  difficulty 
of  the  isocline  method  lies  in  the  labor  required  to  get 
the  isoclines.   So  if  we  can  use  the  digital  computer  to 
generate  a  program,  call  upon  it  to  draw  the  isoclines  and 
phase  trajectories  on  the  same  phase-plane,  then  the  iso- 
cline method  is  more  useful  and  general  than  other  methods, 

For  time-varying  systems,  the  Delta  method  and  Pell's 
method  are  useful.   Deeskshatulu" s  transformation  method 
(by  using  a  new  plane)  is  a  useful  supplement. 

97 


III.   THE  APPLICATION  OF  PHASE-PLANE  METHODS 

Phase-plane  analysis  of  nonlinear  system  problems  not 
only  gives  a  convenient  display  for  interpreting  results, 
but  also  permits  use  of  different  techniques  for  shaping 
the  phase  trajectory  to  make  it  satisfy  design  requirements. 
So,  it  can  be  said  that  the  phase  plane  has  become  a  logic- 
al tool  in  the  solution  of  problems. 

Generally,  phase-plane  analysis  is  limited  to  second- 
order  systems.  Since  higher-order  derivatives  cannot  be 
plotted  on  the  phase  plane,  higher-order  systems  would  not 
be  completely  defined.  Of  course,  the  method  might  be  ex- 
tended from  phase  plane  to  phase  spaces  of  n  dimensions  to 
handle  nth-order  systems,  but  then  it  loses  some  of  its 
simplicity  and  convenience. 

In  this  chapter,  the  manipulations  of  the  phase  plane 
are  discussed,  then  the  uses  of  different  graphical  methods 
and  techniques  for  phase-plane  analysis  of  nonlinear  system 
problems.   All  the  problems  are  of  second  order  or  can  be 
approximated  by  second-order  differential  equations. 

A.    MANIPULATIONS  ON  THE  PHASE  PLANE 

1.   Elementary  Algebraic  Manipulations 
a.   Addition 

Addition  of  phase-plane  trajectories  is  accom- 
plished in  a  fairly  straightforward  manner.   Suppose  we  have 
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two  quantities  x.  and  x„,  and  it  is  desired  to  find  the 
phase-plane  trajectory  of  x,  +  x~. 

Let  x  and  x  represent  phase-plane  coordinates: 

therefore 


Xt  "  Xl  +  X2 


Xt  ~  Xl  +  X2 


(4-1) 


On  the  graph,  simply  add  the  x's  to  obtain  the  new  x  at  a 
particular  time  and  the  x*s  to  obtain  the  new  x.   The  time 
for  each  plot  must  be  known  at  each  point  so  that  addition 
takes  place  at  points  of  isochronism. 
b.   Multiplication: 


If  the  trajectories  are  x,  and  x„  then  the 


products  are: 


Xt  =  Xl  *  X2  (4-3) 


x.  =  X-.X,,  +  XpX-  (4-4) 


Again  the  addition  and  multiplication  must  take  place 
at  points  of  isochronisnu 

c.   Raising  to  a  power: 
For  nth  power 

xfc  =  xn  (4-5) 

x.  =  nx   x  (4-6) 


The  plot  is  taken  point  by  point  and  the  x  coordinate 
is  raised  to  the  nth  power  to  obtain  the  new  x  ;   the  x  co- 
ordinate is  raised  to  the  x    power  and  multiplied  by  nx 

in  order  to  obtain  the  new  x,  . 

t 

d.   Operation  on  the  trajectory  with  a  function: 
If  the  trajectory  for  x  is  known,  then  the 
trajectory  for  f(x)  can  be  founded  as  follows: 

xfc  =  f(x)  (4-7) 


_  df (x)   dx 
Xt  :   dx     dt 


df  (x) 
dx 


(4-8) 


The  new  x  is  the  value  of  f(x)  and  the  new  x  is  the 
slope  of  the  f (x)-vs-x  curve  times  x. 

2.   Differentiation 

The  trajectory  x  is  known  and  the  trajectory  of 
the  derivative  of  x  is  desired.   This  has  a  practical  appli- 
cation.  Suppose  a  device  or  system  output  is  to  be  display- 
ed as  a  phase-plane  trajectory.   This  can  be  done  by  differ- 
entiating  the  output  x,  thus  obtaining  -rr     or  x.   Then  x 
is  plotted  against  x.   Unfortunately,  it  is  very  difficult 
to  differentiate  properly;  and  alternate  is  to  take  the  in- 
tegral of  the  trajectory,  then  differentiate  the  trajectory 
graphically  to  obtain  the  true  x-x  plot. 
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For  illustration  purposes,  let 

u  =  fxdt  ;  (4-9) 


where 


Then 


But 


du 
x  "  dt 


*  _  dx 
x  ~  dt 


dx  . .  dx   du  .  dx 
dt    du   dt    du 


So  in  order  to  find  the  new  x  for  each  point  x  and  u: 

dx 

a.  Find  the  slope  -r—  . 

b.  Multiply  this,  value  by  the  value  of  x  at  that 
point. 

c.  Plot  the  new  x  on  the  -u  axis.  (See  Fig.  3-1) 
3.   Transformation 

It  is  possible  to  transform  the  phase-plane  tra- 
jectory by  applying  some  transformation  to  x  and  x  quanti- 
ties for  some  purpose. 
For  examples 

Xp    r  vX-j  g    x.  j 

X2  =  g^Xl'  Xl^ 
An  example  of  this  might  be 

x~  =  x.  sin  x, 

X~  =  x,  cos  x, 
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-V 


Fig.  3-1.   Differentiation  procedure 
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4.   Methods  of  Shaping  Phase-Plane  Response 

There  are  several  ways  in  which  a  phase-plane  re- 
sponse can  be  shaped  to  conform  to  design  requirements. 
These  are: 

(1)  Choosing  specific  nonlinearities. 

(2)  Limiting. 

(3)  Piecewise  shaping. 

a.  Choosing  Specific  Nonlinearities: 

For  the  following  second-order  differential 
equation: 

x  +  f(x,  x)x  +  g(x,  x)  =  0 

various  combinations  f(x,x)  and  g(x,x),  may  be  chosen  to 
obtain  a  required  response. 

b.  Limiting: 

The  values  of  x  and  x  may  be  restricted  to 
certain  values  as  shown  in  Fig.  3-2.   The  limiting  may  be 
extended  to  x  also  or  even  to  some  function  of  x  and  x. 

c.  Piecewise  Shaping: 
See  Fig.  3-3. 

B.    COMMON  PHYSICAL  NONLINEARITIES 

1.   Saturation:  (the  most  common  of  all  nonlinear 
phenomena) 

Consider  the  case  of  a  saturating  electronic  ampli- 
fier driving  a  motor  in  an  instrument  servo.   (See  Fig. 
3-4 (a)  and  (b) ) 
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Fig.    3-2.      Limiting   of   x   on    the   phase-plane, 


*-/  f(x)+  X 
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/+  £6t'/  +  x=o 
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X+fy*)+K-0 


*   +f5 


(*)  4-  X'O 


x  +  f6(M+x-o 


Fig.  3-3.   Piecewise  shaping  of  trajectory, 
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Amplifier 

]';';r  and  gears 
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Linear  gain  =  X 

(a)   Block  diagram  of  a  servo 
with  saturation. 


Output . 

Linear        ^ 

approximation   ^ 


►  Input 


(b)   Input  vs  output  characteristics. 


Fig.  3-4.   Common  physical  nonlineari ties  -  Saturation, 
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The  equation  in  the  linear  region  is 


i      K  K  K 

E  +  -£-   E  +  -^ — -  E  =  0 
T        NT 
m         m 


-E  <  E  <  +  E 


For  saturation  operation: 


1  • 
E+^^E=+c        E<-En 

m 


E+-^E  =  -c        E^+En 
T  1 

m 


where 


KK  K 
-    s  m  _ 

C  "   NT    *!  ' 

m 


Using  the  isocline  method,  draw  isoclines  on  the  phase- 
plane  E  vs  E.   In  the  linear  region  the  isoclines  are  radial 
straight  lines  (see  Fig.3-4c),  emanating  from  a  focus  at 
the  origin. 

While  in  the  saturation  region  the  isoclines  are  hori- 
zontal, the  transition  from  linear  to  nonlinear  operation  is 
indicated  by  a  vertical  line  at  E  =  E..   Dividing  lines  on 
the  phase  plane  may  be  obtained  from  many  nonlinear  phenom- 
ena. 
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Fig.  3-4. 


(c)   Isoclines,  dividing  lines  and 
trajectory. 
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Fig.  3-5.   Dead  zone  -  (a)  Input  vs  output  character- 
istics (b)  Block  diagram  of  a  servo  with 
dead  zone  (c)  Isoclines,  dividing  lines 
and  trajectory. 
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2 .   Dead  Zone 

This  arises  in  mechanisms  which  are  spring-loaded 
to  minimize  backlash  and  in  many  other  devices  which  are 
insensitive  to  small  signals.   Consider  a  servo  with  a  dead 
zone  in  the  error  transducer. 


E+  —  E  +  E  =  0 
Tm 


-  E..  <  E  <  +  E 


.       K  KK 

E  +  —  E  +  -| — -(E  -  E..)  =  0 
T         Nt  1 

m        m 


E  |  >  E, 


In  the  dead  zone,  there  is  no  control  forcing- function, 
so  in  this  case  isoclines  are  straight  lines  with  slope 

•  1 

-E  -E  f- 

M  =  53- 


(See  Fig.  3-5) . 

3.   Backlash 

Most  mechanical  linkages  exhibit  a  certain  amount 
of  backlash.   (See  Fig.  3-6) 

Consider  a  simple  positioning  servo  with  backlash. 
When  the  backlash  is  taken  up,  the  operation  is  linear.   The 
equation  of  the  system  is 

E  +  2tiX>   E  +  Ui   2E  -  0 
n      n 
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i .  Output 
position 


(a) 


(b) 


M-0,25    M =0-5 


=-0.2 


(c) 


Fig. 3-6.   Common  physical  nonlineari ties 
Backlash. 
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When  zero  velocity  is  reached  and  the  motor  reverses,  the 
output  shaft  remains  stationary  until  the  backlash  is  taken 
up.   Thus  the  error  remains  fixed  at  some  value  E. . 
The  equation  of  this  motion  is 

B     +   2Tu   0  =  En 

m      n  m     1 

4.   Instrument  Servos  with  Coulomb  Friction  and 

Stiction 

a.   A  second-order  servo  with  coulomb  friction, 
(see  Fig.  3-7 (a) ) . 

The  error  equation  for  a  step-displacement  input  is 

j£  +  fE  +  KE  +  C  Sgn(E)  =  0 

from  which  the  isocline  equation  becomes 


•       KE     C  Sgn(E) 
MJ+f    MJ+f 


Note:   Let   E  =  y 


E  =  *L  =    *£.      dE  =  ME 

15   dt   dE   dt   u 

The  second  term  introduces  the  effect  of  the  nonlinear- 

ity.   This  effect  is  seen  to  be  translation  of  the  focal 

C  C 

point  to   E  =  +  —  when  E  is  negative  and  to   E  =  -  —  when 

.  • 

E  is  positive.   E  =  0  (or  E  axis)  is  a  dividing  line.  (See 
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Fig.  3-7).   The  trajectory  terminates  at  D  because  the 
drive  torque  produced  by  the  error  at  D  is  less  than  the 
coulomb  torque  C. 

b.  If  stiction  is  present  the  isocline  equation 
becomes 

•  _     KE     S  Sgn(E) 
MJ+f    MJ+f 

s 

This. moves  the  focal  point  to  -  —  .  (See  Fig.  3-7 (b) ) . 

is. 

A  small  initial  disturbance  (such  as  might  be  represented 
at  point  A)  would  not  produce  sufficient  torque  to  pull  the 
system  free  from  the  stiction.   For  a  large  disturbance 
(such  as  A")  the  system  does  pull  free  and  begins  to 
acquire  a  velocity. 

c.  When  a  ramp  input  is  applied  and  coulomb  fric- 
tion is  present  but  stiction  may  be  neglected,  the  isocline 
equation  becomes 

fOJ.  +  C  Sgn(0   -  E) 
E  =  -  r^-r  + 


MJ+f         MJ+f 

where  U).    =    8   is  the  input  ramp  velocity  and  the  coulomb 
friction  is  designated  by  C  Sgn(0R-E)  =  C  Sgn ( 0  )  because 
the  coulomb  friction  effect  is  controlled  by  the  relative 
velocity  between  the  friction  surfaces  at  the  output. 

The  focus  is  plit  by  the  C  Sgn(0  -E)  term,  and  is 
translated  by  both  C  and  f  U).  .   The  foci  thus  appear  at 

E  =  (fa*.  ^  C)/K;   they  are  displaced  i  C/K  from  the  point 

fO). 
l 

K   * 
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Fig.  3-7.   Phase  trajectory  for  an  instrument 
servo  with  Coulomb  friction  and 
stiction:   Step  response. 
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The  dividing  line  is  a  horizontal  straight  line  at 


Note  j 


at 


E  =  oj. 


Qr>    ~   E  =  8^ 
R         C 


ec  =  o  ,    ft  =  eR  =  « 


(See  Fig.  3-8  (a)).   For  a   ramp  input  suddenly  applied  at 
t  =  0,  the  initial  conditions  are  E(0)=0,E(0)  =  03.,  which 
correspond  to  point  P  on  the  figure. 

d.   When  stiction  is  present  but  discontinuous, 
the  isocline  system  of  Fig.  3-8 (a)  is  not  altered.   But  the 

phase  trajectory  cannot  break  free  from  the  dividing  line 

C  S 

at  E  =  —  .   It  must  remain  on  the  dividing  line  until  E  =  — . 

(See  Fig.  3-8 (b) ) . 

C.    APPLICATIONS  OF  THE  ISOCLINE  METHOD 

Since  there  are  many  methods  of  solution  of  differen- 
tial equations  which  are  comparatively  simple  on  the  phase- 
plane  but  difficult  or  impossible  otherwise,  it  is  rather 
obvious  that  the  phase-plane  provides  a  very  convenient 
means  of  analysis  of  extremely  nonlinear  first  and  second- 
order  differential  equations.   It  seems  that  the  isocline 
method  is  a  general  and  useful  method  with  which  to  approach 
problems  by  using  the  phase-plane  method.   Application  of 
the  isocline  method  to  specific  problems  illustrates  the 
usefulness  of  this  method. 
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Fig.  3-8a. 


Phase  portrait  of  an 
instrument  servo  with 
Coulomb  friction  in 
response  to  a  ramp 
input. 


Fig.  3 -8b. 


Effect  of  stiction  on 
the  ramp  response  of  a 
servo  with  Coulomb 
friction. 
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1.   First-Order  System 

The  first-order  equation  has  the  form 

ff  =  f(x,t)  (3-1) 

Assume  that  the  function  f(x,t)  is  continuous  and  single- 
valued  with  the  possible  exception  of  certain  singular 
points.   The  isocline  method  requires  a  graphical  construc- 
tion performed  on  axes  of  x  and  t  and  yields  a  solution  as 

a  curve.   For  any  given  point  on  the  x-t  plane,  the  numeric- 

dx 
al  value  of  f(x,t),  and  thus  of  -rr  ,  can  be  calculated. 

Example  1:   Find  a  solution  for  the  equation 

dx 

—  =  -  x  +  t 

dt 

Consider  the  initial  conditions  x   =0   at   t   =0,  and 

o  o 

also   x   =  0.2   at   t   =0.   For  this  linear  equation,  the 
o  o  •*       » 

isoclines  are  straight  lines  given  by  the  relation  m  =  x  +  t. 
A  family  of  isoclines,  each  carrying  its  line  segments  of 
proper  slope  is  shown  in  Fig.  3-9.   Solution  curves  are 
sketched  in  the  figure,  starting  at  each  of  the  specified 
initial  conditions. 

Example  2:   RL  circuit  with  nonlinear  L 

/      R 

S    /yV\A. 


—  E 


L  ~  +  iR  =  Eu(t) 
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Since  the  inductor  is  nonlinear,  and  the  coefficient  L 
therefore  variable,  it  is  more  convenient  to  express  the 
relationship  as 


M^+    iR   =    Eu(t) 
dt 


from  which 


a<b   =   Eu(t)     -    iR 
dt  M 


The  equation  of  the  isoclines  is  obtained  by  letting 

-r-r   =  M,  some  constant,  and  noting  that  Eu(t)  =  E   for  t  >  0, 

then 

E  -  MM., 
i  =  ^ 


Knowing  the  magnetization  curve  of  the  inductor,  one 
can  obtain  the  transient  current  variation  of  the  nonlinear 
RL  circuit,  i  vs  t.  (See  Fig.  3-10) 

2.   Second-order  System 

(Example  1)   The  equation  for  a  damped  pendulum  is 

d  x  ,  T.  dx  ,   .       _ 

+  K  -tt   +  sin  x  =  0 


Let 


dt       dt 


dx 
v  =  dt 


2 

d  x    dv    dv  m    dx  _  dv 

,,2    dt    dx    dt    dx 
dt 


v 
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0.1         0.2        0.3        0.4         0.5        0.6         0.7        0.8 

Fig.    3-9.        (Example    III-C-1-A) . 
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Fig.    3-10.       Example    (III-C-1-B) . 
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Thus 


d  x  ,  T/  dx  . 

7T  +   K  -rr  +  sin  x 

,,2      dt 
dt 


-=—  +  Kv  +  sin  x  =  0 
dx 


or 


dv      (Kv  +  sin  x) 
dx  v 


Let 

K  =  1 


dv  _    (v  +  sin  x) 
dx         v 


j—  is  the  slope  at  any  point  x,  v  of  a  possible  phase- 
plane  trajectory. 

Now  assume  a  number  of  values  for  -r— 

dx 

-t—  Equation  of  v  vs  x 


0  v  =  -  sin  x 

1  v  =  -  y  sin  x 
-1  0  =  -  sin  x 

2  v  =  -  -r  sin  x 

-2  v  =  sin  x 

1   ■ 

3  v  =  -  t  sin  x 


4 
1 

2 


-3  v  =   77  sin  x 


(See  Fig.  3-11) 
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Fig.    3-11.       Isocline    solution. 
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Example  2%      Nonlinear  system  x  +  xx  +  x  =  0 

Let   x  =  Mx  ; 

then 

Mx  +  xx  +  x  =  0 


-x 

x  = 


M  +  x 

Draw  a  number  of  straight  lines  x  =  Mx  through  the 
origin  (corresponding  to  different  values  of  M)  on  the 
x-x  plane  (as  in  Fig.  3-12 (a)). 

Plot  x  vs   x  from  equation  x  +  xx  +  x  =  0  for  con- 
stant (different)  values  of  x  (see  Fig.  3-12 (a)). 

With  any  particular  straight  line  in  view,  say  WL  —   M_, 
determine  all  possible  points  of  intersections  of  above 
plots.   The  procedure  is  repeated  for  M  =  M?#  M_,  etc., 
until  a  fair  number  of  isocline  curves  are  obtained  on  the 
phase  plane.   The  respective  slope  lines  are  marked  on 
these  isoclines.   Joining  these  slope  lines  gives  the  dif- 
ferent trajectories  or  the  entire  phase  portrait  (as  in 
Fig.  3-12 (b)). 

Example  3:  Consider  the  case  of  a  saturating  elec- 
tronic amplifier  driving  a  motor  in  an  instrument  servo. 
(See  Fig.  3-13 (a) ) 

The  equation  in  the  linear  region  is 

,      K  KK 

E  +  ~  E  +  -~ E  =  0       -E.<  E  <  +E.. 

T        NT  11 

m        m 
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t^'-to 


*=-.* 


Fig.    3-12 (a)      Example   x   +  xx   +  x   =   0 
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Fig.    3-12 (b) .      Example      x   +  xx   +  x   =   0 
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Fig.    3-13 (a) 
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where  E   is  the  value  of  |  E  |  which  causes  saturation  and  N 
is  the  gear  ratio. 

For  saturation  operation,  the  equations  are; 


E  +  it-   E  =  +  C 
T 
m 


E  <  -E. 


E  +  ^-   E  =  -  C 
T 
m 


E  £  +E. 


where 


C  = 


K  K  K 
s   m 


NT 


m 


On  the  E  vs  E  phase  plane  the  isoclines  are  all 
straight  lines;   in  the  linear  region  the  isoclines  are 
radial  straight  lines  emanating  from  a  focus  at  the  origin, 
While  in  the  saturation  region  the  isoclines  are  horizon- 
tal, the  phase  trajectory  is  constructed  in  the  usual 
manner  (as  in  Fig.  3-13 (b) ) . 

Example  4s   Oscillator  with  relay  control. 


/  /  /  // 


mx  +  Kx  =  F 


K      F 

x  +  —  x  =  — 

m     m 


Let 


£-  i 

m 


F  A 

—  =  u 
m 


then  x  +  x  =  u 
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Let 

x  =  y 


then 


where 


dt    dx    dt    Mx 


dx 


thus 

My  +  x  =  u 


u  -  x 

M  =  — 


a.  If  u  =  K  the  phase  trajectory  is  as  in 
Fig.  3-14 (a). 

b.  If  u  =  -K  Sgn(x) 


M  =  - 


x  f  K  Sgn(x) 


x  =  0    or   y-axis  is  a  dividing  line. 
(See  Fig.  3-14 (b) ) 

x  "I"  K 

If  x  is  positive   M  =  - 


if  x  is  negative   M  =  - 


y 

x  -  K 


c.   If  u  =  -K  Sgn(y) 

M  =  _  x  +  K  Sgn(y) 

y 

y  =  0   or   x-axis  is  a  dividing  line, 
(See  Fig.  3-14(c) ) 
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Lr 


-9 


-15 


Fig.  3-13 (b) 


15 

9 

5 
4 


X-Scale 
Y-Scale 


2.00E-01  Units  Inch. 
2.00E-01  Units  Inch. 
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If  y   is  positive   M  =  - 
If  y  is  negative    M  =  - 


X     + 

K 

Y 
x    - 

K 

Y 


c.   If  u  =  -K  Sgn(x+y) 


M  = 


x  -  K  Sgn(x+y) 
Y 


x  +  y  =  0   is  a  dividing  line, 
(See  Fig.  3-14 (d) ) 

e.   If  u  =  -K  Sgn(x-y) 


=  x  -  K  Sgn(x-y) 

Y 


x  -  y  =  0   is  a  dividing  line. 
(See  Fig.  3-14(e) ) 

Due  to  the  changing  dividing  line  the  system  became 
unstable. 

Example  5 : 


u 


(s+2) (s+  |) 


x 


Given: 


u   =' 


x   <  0.5 


-Sgn(x)  ,    I  x  I  ^  0.5 
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Since 


x  (sj 
u(s) 


(s+2)  (s+  -|)    s2  +  |-  s  +  1 


Let 


then 


5  • 
x+-jx+x  =  u 


x  =  y 


dy    dy    dx 
X  =  -rr  =  -T-    "    -rr   =    MX 
dt    dx    dt 


thus 


M  = 


-  2 


y  -  x  +  u 


At   u  =  0,  the  equation  of  isoclines  is 


My  =  -  2   Y 


-  t  y  -  x 


At 


u  =  -  Sgn(x) 


If  x  is  positive,  the  equation  of  isoclines  is 


My  =  -  -J  y  -  x  -  1 


If  x  is  negative,  the  equation  of  isoclines  is 


My=--~y-x  +  l 


(See  Fig.  3-15) . 


128 


X-Scale  =  l.OOE+00  Units  Inch. 
Y-Scale  =  l.OOE+00  Units  Inch. 


Fig.  3-14(a) 
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■J-Im*;  ,t  ■ 


X-Scale  =  2.00E-01  Units  Inch. 
Y-Scale  =  2.00E-01  Units  Inch. 


Fig.  3-14 (b) 
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/ 


Y 


Isocline  is  the  sam 
as  3-14 (b) . 


\ 


-•• 


■»■*• 


/  dividing  line 


*•>* 


Fig.  3-14 (c) . 
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Isocline  is 
the  same  as 
3-14(b) . 


Fig.  3-14 (d) 
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I,  Y 


Isocline  is 
the  same  as 
3-14(b) . 


Fig.    3-14 (e) 


133 


Fig.  3-15 
X-Scale  =  5.00E-01  Units  Inch. 
Y-Scale  =  5.00E-01  Units  Inch. 
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3.   High -Order  Systems 

The  phase-plane,  as  normally  used,  is  restricted 
to  problems  which  can  be  described  by  first  and  second- 
order  differential  equations.   This  is  a  fairly  serious  re- 
striction, since  many  practical  systems  require  higher-order 
differential  equations  for  an  adequate  description.   In  the 
discussion  of  numerical  methods  for  solving  differential 
equations,  first-order  equations  are  studied  in  the  begin- 
ning.  Then  equations  of  higher  order  are  solved  by  replac- 
ing each  equation  of  nth  order  by  n  equations  of  first-order 
and  solving  these  simultaneously.   In  like  manner,  the  dis- 
cussion of  graphical  methods  has  begun  by  considering  the 
isocline  method  of  solving  a  first-order  equation  first, 
then  second-order  equations.   But  further  consideration  may 
not  be  possible  since  there  is  not  enough  information  avail- 
able initially  to  construct  isocline  curves  for  the  several 
first-order  equations  that  would  arise  from  a  single  high- 
order  equation.   So,  if  we  can  factor  high-order  equations 
into  several  second-order  equations,  these  can  be  solved  by 
phase-plane  method.   The  difficulty  is  in  converting  a 
high-order  equation  into  several  second-order  equations. 

No  convenient  long  hand  method  is  available,  but  use 
of  the  digital  computer  permits  study  of  such  phase  space 
problems. 
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IV.   GENERATING  THE  PROGRAMS  TO  DRAW  ISOCLINES  AND 
TRAJECTORIES  OF  SECOND-ORDER  DIFFERENTIAL 

EQUATIONS 

A.    x  +  F(x)  x  +  G(x)  =  0  type  equations  (See  Appendix, 
Program  1)   (Where  F(x)  f  G(x)  can  be  any  function  of  x) 

x  +  F(x)  x  +  G(x)  =0 


Let 


Let 


thus 


y  =  x 


dy  .  dy_   dx 
dt    dx    dt 


¥   =  M 

dx 

then 

x  =  My 

x  +  F(x)  x  +  G(x) 
=  My  +  F(x)  y  +  G(x)  =0 

v  =  -  — ^B — r  equation  of  isoclines 

1  M  +  F(x)    ^ 

Trajectories : 

ZD0T(1)  =   Z(2) 

ZD0T(2)  =  -F(Z(1))*  Z(2)  -  G(Z(1)) 

(By  using  subroutine  program  RKLDEQ.) 
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Example    1:      x   +  xx   +   x   =   0 


li*uti 


Fig.  4-1.     X-Scale  =  l.OOE+00  Units  Inch. 
Y-Scale  =  l.OOE+00  Units  Inch. 
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*2 

Example  2:    x  -  e(l-x  )x  +  x  =  0 


Ml 


Fig.  4-2.     X-Scale  =  8.00E-01  Units  Inch 
Y-Scale  =  8.00E-01  Units  Inch 
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Example    3:        x"   +(cos   x)x   +    tan  x 


rH 


0 


Fig.    4-3.  X-Scale   =    l.OOE+00   Units    Inch, 

Y-Scale   =    l.OOE+00    Units    Inch, 
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Example  4:    x+25(l+0.1x)=0 


Fig.  4-4.     X-Scale  =  l.OOE+00  Units  Inch, 
Y-Scale  =  l.OOE+00  Units  Inch, 
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Example  5:    x  +  0.5x  +  x  =  0 


-3 


2b 


Fig.  4-5.     X-Scale  =  l.OOE+01  Units  Inch. 
Y-Scale  =  l.OOE+01  Units  Inch. 
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Example  6:    x  -(l-x2)x  +  x  =  0 


Fig.  4-6.     X-Scale  =  l.OOE+00  Units  Inch. 
Y-Scale  =  l.OOE+00  Units  Inch. 
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Example  7  :    x  +  x  +  x  |  x  |  =0 


Fig.  4-7.     X-Scale  =  2.00E+00  Units  Inch. 
Y-Scale  =  2.00E+00  Units  Inch. 
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Example  8:    x+  —  x  +  x  =  0 


Fig.  4-8.     X-Scale  =  l.OOE+00  Units  Inch. 
Y-Scale  =  l.OOE+00  Units  Inch, 


144 


Example  9:    x  +  x  +  sin  x  =  0 


i 


Fig.  4-9 


X-Scale  =  9.00E-01  Units  Inch, 
Y-Scale  =  9.00E-01  Units  Inch, 
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Example    10:        x   +    (1-  ]  x  |  ) x   +  x   =   0 


99 

?9 

59 
39 


:  n 


ifl 


Fig.    4-10. 


X-Scale  =  1.00E+00  Units  Inch, 
Y-Scale  =  1.00E+00  Units  Inch, 
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Example  11 


x  + 


sin  x  =  0 


50 


59 

39 


Fig.  4-11 


X-Scale  =  2.00E+00  Units  Inch. 
Y-Scale  =  2.00E+00  Units  Inch. 
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2 
Example  12:    x  +(l-x  )x  +  x  =  0 


Fig.  4-12 


X-Scale  =  5.99E-01  Units  Inch. 
Y-Scale  =  5.00E-01  Units  Inch, 
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Example  13:    x  +  2  xx  +  x  =  0 


g 


Fig.  4-13. 


X-Scale  =  l.OOE+00  Units  Inch. 
Y-Scale  =  l.OOE+00  Units  Inch. 
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•  2 

B.    x  +  x   +G(x)=0  TYPE  EQUATIONS  (SEE  PROGRAM  2) 

(WHERE  G(x)  IS  ANY  FUNCTION  OF  x) 

x  +  x2  +  G(x)  =0 


Let 


then 


thus 


x  =  y 


dy   dx   .. 

x  =  -&     dt  =  My 


x  +  x   +  G(x) 


=  My  +  y2  +  G(x)  =  0 


y2  +  My+  G(x)  =  0 


-M  +  j¥?-   4G(x)  ^  ^ 

Y  2  ^ 

— 

a  j  ^  ^equations  of 

-"    isoclines 


-MZ 


-M  -  M    -    4G(x) 
Y  =  5 


Trajectories : 

ZD0T(1)  =  Z(2) 

ZD<£T(2)  =-Z(2)  **  2  -  G(Z(1)) 
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Example  14:    x  +  x   +  x  =  0 


-3 


Fig.  4-14 


X-Scale  =  5.00E-01  Units  Inch, 
Y-Scale  =  5.00E-01  Units  Inch, 


151 


•  2  2 

Example    15:         x+x      +x      =0 


Fig.  4-15 


X-Scale  =  5.00E-01  Units  Inch, 
Y-Scale  =  5.00E-01  Units  Inch, 
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.2 


Example    16:         x   +  x      +x|x!     =0 


Fig.  16. 


X-Scale  =  5.00E-01  Units  Inch, 
Y-Scale  =  5.00E-01  Units  Inch, 
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C.         |x|x+x   +   x   =   0 

TYPE    OF   EQUATIONS 
x   +   x|x|+x   =   0 

(SEE    PROGRAM   3) 

x    =   x|x1+x    =    0 


In   Section  A   and   Section   B    the   equations   of   isoclines 
are    solved  by  y   in    terms    of  x.      In    this    section  we  will 
solve   x   in    terms   of   y    to  prevent   imaginary  values. 

1.       |  x   |  x   +  .x   +  x   =   0 

Let 

y  =  x   , 

then 

*  =  1^  It  =  M^ 

thus 

|x   |x   +  x   +  x   =  |y|My   +  y   +  x   =   0 

x   =   -  |  y  |  My   -   y 

=   -  y(   I  y  1  M  +   1)  equation   of   isocline. 

Trajectories : 

ZD<£T(1)    =   Z(2) 

zDrtTm  =    "z(2)  -  z(1> 

ZD0TU)     -         DABS     (Z(2)) 
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2.      x   +x|x|+x   =   0 

x+x|x|+x=My+x|y|+x=0 
My 


x  =  - 


equation  of  isoclines 


1    +   I  y   I 

Trajectories : 

ZD0T(1)    =   Z(2) 

ZD0T(2)     =-Z(l)*    DABS(Z(2))     -    Z(l) 

3.       x    +   x    |  x    |     +x    =   0 

x   +  x  |  x  |     +x=My+y|y|     +x=0 

x   =    -  My   -   y  1  y  1  =   -  y (M  +  1  y  |  )       equation   of 

isoclines. 

Trajectories : 

ZD0T(1)    =   Z(2) 

ZD0T(2)    =-Z(2)*    DABS(Z(2))     -    Z(l). 
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Example    17:         x   +  x|x|+x   =   0 


Fig.  4-17. 


X-Scale  =  l.OOE+00  Units  Inch, 
Y-Scale  =  l.OOE+00  Units  Inch, 
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Example    18:         x   +   xlx|+x   =   0 


Fig.  4-18, 


X-Scale  =  l.OOE+00  Units  Inch, 
Y-Scale  =  l.OOE+00  Units  Inch, 
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Example    19:         |x|x   +   x+x   =   0 


Fig.  4-19 


X-Scale  =  l.OOE+00  Units  Inch, 
Y-Scale  =  l.OOE+00  Units  Inch, 
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D.    USE  OF  PROGRAMS 
1.   Table  1. 


Type  of  Equation 

How 

to  Set  the  Subroutine  and 

Function  Programs 

x  +  F(x)x  +  G(x)  =0 

(1) 

Need  to  compare  standard  form 

(Refer  to  Section  A) 

to  give  real  function  F(x)  and 
G(x)  in  each  function  subpro- 
gram. 

(2) 

Rewrite  the  given -equation 
x  +  F(x)x  +  G(x)  =  0  into  a 
system  of  equations  in  sub- 
routine program. 

ZD0T(1)  =  Z(2) 

ZD0T(2)  =-F(x)*Z(2)-  G(x) 

(F(x),  G(x)  should  change  to 
z  form) 

2 
x  +  x   +  G(x)  =  0 

(1) 

Need  to  give  real  function 

(Refer  to  Section  B) 

G(x)  in  function  subprogram. 

(2) 

Same  as  above  (2) . 

1  x  1  x  +  x  =  0 

(1) 

Isoclines : 

x  +  x|x|+x  =  0 

• 

Let   x  =  y 

x  +  x|x|+x  =  0 

solve x in  terms  of  y  and  give 

(Refer  to  Section  C) 

this  equation  in  main  program. 

(2) 

Same  as  above  (2) . 

Table  .1. 
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Input  deck,  See  Table  2. 


Sequence  of  Cards 

Purpose 

Punch  in  Column 

1st  Card 

Scale 

( Number  of  Cards 
<of  different 
granges  of  slope 

One  significant 
number 

Column  1-10 
(column  11-20 

2nd  Card 

Title 

Column  1-48 

3rd  Card 

Title 

Column  1-48 

4th  Card 

• 
• 

N 
No.  of  cards 
depends  on 
number  of  dif- 
ferent ranges 
of  slope  cards. 

Different  ranges 

Column  1-10  Ini- 
tial 
value 

Column  11-20 

Increment 

Column  21-30  Final 
value 

(N+l) th  Card 

Time  increment 

Number  of    ) 
trajectories  j 

Column  1-10 
Column  11-20 

(N+2) th  Card 

• 
• 

K 

No.  of  cards 
depends  on 
number  of  tra- 
jectories. 

Initial  time 
Final  time 
Initial  x 
Initial  y 

Column   1-10 
Column  11-20 
Column  21-30 
Column  31-40 

(K+l) th  Card        Blank 

! 1 

Blank 

Table  2. 
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NOTE:   FOR  SHIFTING  THE  CENTER  OF  COORDINATES.  SEE 


AND  COMPARE  PROGRAM  10. 

4-  RR 
Example  20.   x  +  J^63  *  cos  x  +  sin  x  =  ° 


!    1 


i      hH-i  -    - 


l. 


w 


■\\ 


n. 


!      i 


\\Y\\    Uito 


^pfcg»  ■„« 


Fig.    4-20. 


X-Scale  =  2.00E+00  Units  Inch. 
Y-Scale  =  2.00E+00  Units  Inch. 
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E.    DISCUSSION 

In  Chapter  II  it  was  shown  that  the  isocline  method  is 
a  very  useful  general  method  of  phase-plane  analysis,  the 
only  difficulty  being  the  labor  required  to  get  the  iso- 
clines.  This  difficulty  has  been  reduced  by  developing 
computer  programs  to  draw  both  isoclines  and  trajectories. 

In  Section  A,    a  program  works  for  x  +  F(x)x  +  G(x)=  0 
type  equations.   This  is  the  general  case,  and  also  can  be 

done  by  hand  labor.   In  Section  B,  a  program  works  for 

•  2 
x  +  x   +  G(x)  =  0  type  equations.   Working  this  kind  of 

equations  by  hand  labor,  some  supplementary  methods  are 
needed,  such  as  Deekshatulu ' s  or  Murthy's,  but  use  of  the 
computer  saves  time  and  labor. 

In  Section  C,  a  program  works  for  some  special  prob- 
lems.  In  this  program  the  equation  of  isoclines  is  first 
solved  for  x  in  terms  of  y  (Note:   In  Section  A   and  B,  iso- 
cline equations  are  solved  for  y  in  terms  of  x) .   After 
calculating  the  required  points  the  plotting  routine  treats 
y  as  a  function  of  x„   The  reason  for  this  is  to  avoid  the 
problem  of  complex  values  in  the  computation  of  y.   This 
procedure  is  not  always  successful;  there  are  at  least  two 
problems  for  which  complex  values  apparently  cannot  be 
avoided.   They  are: 

1.  x  +  x  I  x  I  +x   =0 

2.  x  +  xlx    +  x  |x  |  =  0 


162 


For  example:   Equation   x  +  x  |  x  |  +  x   =0 


Let 

• 

X 

= 

y 

then 

X 

= 

My 

thus 


x  +  x  I  x  1  +  x   =  My  +  y  I  y  I  +  x   =0 


Case  a:   Solving  for  y  in  terms  of  x 


2         2 

y   +  My  +  x   =0 


Jy? 


-  M  ±  M  4  x 

y  = 


when 

2 

M  <  4x 

y  is  complex. 

Case  b:   Solving  for  x  in  terms  of  y 

x   =  -  [y  |  y  |  +  My] 

x  =  J-[y   I  y  I  +  My] 
when  [y  |  y  |  +  My]  is  positive  x  becomes  complex, 
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V,   GENERATING  PROGRAMS  TO  DRAW  ISOCLINES 
AND  TRAJECTORIES  WITH  DIVIDING  LINES 
OF  THE  SECOND-ORDER  DIFFERENTIAL 
EQUATIONS 

SATURATION 

1.   Program  for  Vertical  Dividing  Lines  only 

(See  Program  4) 

Refer  to  Chapter  III,  Section  B-l. 
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Example  1.    Tra jec tories :   x  +  0 . 2x  +  x  =  0 
Vfor  Inner 
regions  )  Dividing  lines 

x  =  ±  .3 

Isoclines : 

/  „    „.    c  n \  Inner  linear  region 

(  See  Fig.  5-1)  ^ 

x  +  0.2x  +  x  =  0 
My  +  0.2y  +  x  =  0 
-x 


y  = 


Trajectories : 


M  +  0.2 


Saturation  regions: 
x  +  0.2x  ±  0.3  =  0 
My  +  0.2y  i  0.3  =  0 
T  0.3 


Y  = 


M  +  0.2 


Inner  region 

ZD0T(1)  =  Z(2) 

ZD<2T(2)  =  -0.2*  Z(2)  -Z(l) 

Saturation  region 

ZDtfT(l)  =  Z(2) 

ZD$T(2)  =  -0.2*  Z(2)  ±  0.3 
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if 


Fig.  5-1. 


X-Scale  =  2.00E-01  Units  Inch. 
Y-Scale  =  2.00E-01  Units  Inch. 
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Example    2:         x+0.2x+x=0 

Dividing  lines   x  =  +  0.2 


Fig.  5-2. 


X-Scale  =  2.00E-01  Units  Inch 
Y-Scale  =  2.00E-01  Units  Inch, 
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2.   Program  for  Vertical  Dividing  Lines  which  can  be 
Rotated.  (See  Program  5) 
Equation  of  dividing  lines: 

Ax  +  By  =  ±  C 
Isoclines  in  linear  region: 

x  +  F(x)x  +  G(x)  =  0 

assume   F(x)  =  1,   G(x)  =  x 
Let   y  =  x  ,  then  x  =  My 

Y  :  M  +  1 

Isoclines  in  saturation  regions: 
±  C 


Y  = 


M  +  1 


Trajectories : 

Inner        fzD0T(l)  =  Z(2) 

region     |^ZD0T(2)  =  -'1    *  Z(2)  -  Z(l) 

a     l.  i  ZD0T(1)  =  Z(2) 

Saturation   ■ 

region     I  ZD0T(2)  =  -Z(2)  ±  C 
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Example  3:    Dividing  lines:   x  +  0.8y  =  ±  0.3 
Trajectories   :   x  +  0.2x  +  x  =  0 


U 


Fig.  5-3 


X-Scale  =  2.00E-01  Units  Inch, 
Y-Scale  =  2.00E-01  Units  Inch, 
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Example  4:    Dividing  lines:   x  -  0 .  8v  =  i  0.3 
Trajectories   :   x  +  0.2>:  +  x  =  0 


Fig.  5-4. 


X-Scale  =  2.00E-01  Units  Inch. 
Y-Scale  =  2.00E-01  Units  Inch. 
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Example  5:    Dividing  lines:   x  =  ±  0.3 

(Y  =  0) 
Trajectories   :   x  +  0.2x  +  x  =  0 


21 


2':. 
?> 

I.' 

-33 


X-Scale  =  2.00E-01  Units  Inch. 
Y-Scale  =  2.00E-01  Units  Inch. 


Fig.  5-5 
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B.    DEAD  ZONE.    (SEE  PROGRAM  6) 
(Refer  to  Chapter  III 
Section  B-2 . ) 

Example  6:        (See  Fig.  5-6) 

Inner  region: 

Isoclines : 

x+0.2x+0.3=0 
My  +  0.2y  +  0.3  =  0 

-0.3 


M  = 


M  +  0.2 


Trajectories : 

(zD<£T(l)  =  Z(2) 
(zD0T(2)  =-0.2*  Z(2) 

External  region: 
Isoclines : 

Right  and  left 

x+0.2x+x^o.3=0 

Trajectories : 

ZD#T(1)  =  Z(2) 

ZD(JT(2)  =  -  0.2  *  Z(2)  -  Z(l)±  0.3 
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Fig.  5-6. 


X-Scale  =  2.00E-01  Units  Inch. 
Y-Scale  =  2.00E-01  Units  Inch. 
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C.    IDEAL  RELAY.    (SEE  PROGRAM  7,  8,  8A) 

For  x  +  0.2x  +  u  =  0 
u  =  -  K  Sgn(x) 

Switching  line  is    x  =  0 

or  y   axis. 

Isoclines : 

At   +x    x  +  0.2x  -  K  =  0 
At   -x    x  +  0.2x  +  K  =  0 

Trajectories : 

At   +x    x  +  0.2x  -  K  =  0 


(zD0T(l)  =  Z(2) 

(  ZD#T(2)  =  -  0.2  *  Z(2)  +  K 


At   -x    x  +  0.2x  +  K  =  0 


(zD0T(l)  =  Z(2) 

(ZD0T(2)  =  -  0.2  *  Z(2)  -  K 
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Example  It         x+0.2x+u=0,     u=0.3  Sgn  x 
Switching  line:   x  =  0   or   y-axis. 

Isoclines : 

-0  3 
Right  of  y-axis:    y  =  M  +'Q<2 


T  CU         4=  •  0.3 

Left  of   y-axis:    y  =  -  +  Q ^ 2 


Trerjec  tories  : 

Right     (  ZD(4T(1)  =  Z(2) 

ZD*T(2)  =  -  0.2  *  Z(2)  +  0.3 


Left      fzD(JT(l)  =  Z(2) 

\zD<£T(2)  =  -  0.2  *  Z(2)  -  0.3 


(See  Fig.  5-7) 
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Fig.  5-7. 


X-Scale  =  2.00E-01  Units  Inch. 
Y-Scale  =  2.00E-01  Units  Inch, 
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Example  8(a) :    x  +  0.02x  +  u  =  0 

u  =  -  0.3  Sgn(x) 


■I 


Fig.  5-8(a) . 


X-Scale  =  2.00E-01  Units  Inch, 
Y-Scale  -  2.00E-01  Units  Inch, 
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Example  8 (b) : 


x  +  0.02x  +  u  =  0 
u  *  -  0.3  Sgn(8x-y) 
(x  =  y) 
+-  -• 


Fig.  5-8(1?) 


X-Scale  -  2.00E-01  Units  Inch. 
Y-Scale  =  2.00E-01  Units  Inch. 
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Example  8(c) 


x  4-  0.0  2x  '-  u  =  0 
d  =  -  0.3  Sgn(8x+y) 
(x  =  y) 


Fig.  5-8(c) . 


X-Scale 
Y-Scale 


2.00E-01  Units  Inch, 
2.00E-01  Units  Inch, 
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Example  5-8 (d) :    x  +  Q.02x  +  u  =  0 

u  =  -  0.3  Sgn(4x-y),  x  =  y 


vj  1 


Fig.  5-8(d). 


X-Scale  =  2.00E-01  Units  Inch 
Y-Scale  =  2.00E-01  Units  Inch 
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Example  5-8 (e) : 


x  +  0.02x  +  u  =  0 

u  =  -  0.3  Sgn  (4x+y) ,  x  =  y 


lig.  5-8(e) 


X-Scale 
Y-Scale 


2.00E-01  Units  Inch. 
2*.00E-01  Units  Inch, 
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Example  5-8 (f) : 


x  +  0.02x  +  u  =  0 

u  =  -  0.3  Sgn(2x-y),  x  =  y 


inj-l 


Fig.  5-8(f). 


X-Scale  =  2.00E-01  Units  Inch. 
Y-Scale  =  2.00E-01  Units  Inch. 
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Rx ample    5-8 (g) 


x   +   0.02x   +    u   =   0 

-    -   0.3    Sgn(2<N      )  ,    x 


=    y 


Pig.    5-3 (g) 


X-Scale  =  2.00E-01  Units  Inch. 
Y-Scale  =  2.00E-01  Units  Inch. 
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Example  5-8 (h) :    x  +  0.02x  +  u  =  0 

u  -  -  0.3  Sgn(x-y),  x  =  y 
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Fig.  5-8(h) 


X-Scale 
Y-Scale 


2.00E-01  Units  Inch, 
2.00E-01  Units  Inch, 
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Example  5-8 (i) :    x  +  0.02x  +  u  =  0 

u  =  -  0.3  Sgn(x+y),   x  =  y 
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Fig.  5-8(i). 


X-Scale  =  2.00E-01  Units  Inch, 
Y-Scale  =  2.00E-01  Units  Inch. 
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Example  5-8 (j) :     x  +  0.02x  +  u  =  0 

u  -  -  0.3  Sgn(y) ,    x  =  y 
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Fig.  5-8(j). 


X-Scale  =  2.00E-01  Units  Inch, 
Y-Scale  =  2.00E-01  Units  Inch, 
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D.    RELAY  WITH  DEAD  ZONE   (SEE  PROGRAM  9) 

For  x+x  +  x  +  u  =  0 
u  =  -K  Sgn(x) 

Isoclines : 

At  dead  zone  region: 

x   +   x   +   x   =   0 

My   +  y   +  x   =   0 


-x 

Y   = 


M   +    1 


At    right: 

x+x+x-K   =   0 
My   +  y   +   x    -K   =   0 


-x+K 
Y 


M   +    1 
At    left: 

x+x  +  x+K  =  0 
My+y+x+K=0 
-x-K 


y 


M  +  1 


Trajectories : 

At  dead  zone  region: 


ZDtfT(l)  =  Z(2) 
x  +  x  +  x  —   0. 

ZD0T(2)  =-Z(2)  -  Z(l) 


87 


At    right: 

(zDdTCL)    =      Z(2) 
x+x+x-K=0,         I 

lzD*T(2)     =-Z(2)     -    Z(l)+K 

At    left: 

(ZD0T(1)     =      Z(2) 

x+x+x+K=0,         ) 

(zD0T(2)     =  -Z(2)-    Z(l)-    K 
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Example   9 : 


x+0.2x+x+u=0 

u   =   -  0.3   Sgn(x) 
dividing    lines:      x   =    ±  0.2 


Fig.  9. 


X-Scale  =  2.00E-01  Units  Inch 
Y-Scale  =  2.00E-01  Units  Inch 


189 


Example  10 : 


x  +  x  +  u  =  0 

u  =  -  0.3  Sgn(x) 
dividing  lines:  x  =  ±0.2 
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Fig.  10 


X-Scale  =  2.00E-01  Units  Inch. 
Y-Scale  =  2.00E-01  Units  Inch. 
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IV.   CONCLUSIONS 

Chapter  III  summarizes  eleven  useful  methods.   Some  of 
them  draw  direction  fields  to  get  trajectories.   Some  of 
them  are  good  supplements  to  the  isocline  method.   Some  of 
them,  by  using  geometrical  relationships,  draw  the  trajec- 
tories step-by-step.   In  general,  none  of  these  methods  is 
a  perfect  one.   For  specific  problems  one  may  be  better  than 
the  others  in  terms  of  the  labor  required.   But  if  we  want 
to  look  at  the  overall  picture  of  some  nonlinear  system 
characteristics  such  as  stability  limit  cycles,  etc.,  none 
of  them  can  do  this  easily.   All  need  lots  of  time  and  hard 
labor. 

The  isocline  method  has  many  advantages.   Programs 
have  been  developed  to  compute  and  draw  both  isoclines  and 
trajectories.   All  the  programs  in  Chapter  IV  and  V,  first 
draw  isoclines,  then  draw  trajectories  from  given  initial 
conditions.   In  analysis  and  design  work,  isoclines  may  not 
be  needed  for  all  problems.   Conversely,  sometimes  only 
isoclines  are  needed  and  we  do  not  need  trajectories.   Use 
of  isoclines  with  superposed  trajectories  gives  a  global 
pitcure.   That  is,  it  indicates  what  can  happen  to  trajec- 
tories for  all  initial  conditions,  and  tells  whether  the 
system  is  stable  or  not,  has  limit  cycle  or  not,  ...  etc. 
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The  programs  generated  in  Chapter  V  add  dividing  lines 
to  the  isoclines  and  trajectories.   By  changing  or  shifting 
dividing  line  position,  the  trajectories  are  shaped.   This 
is  helpful  for  analysis  and  design  work. 

In  this  thesis,  all  graphs  are  plotted  by  machine  and 
all  are  two  dimensional.   The  variable  symbols  y  vs  x  or 
x  vs  x  are  used.   There  are  no  time  variables  in  the  plot. 
If  the  system  equation  involves  time  as  a  variable,  we  may 
choose  the  ft   method  to  work  it  or  apply  some  transforma- 
tion, changing  the  equation  into  the  standard  type  forms. 
In  general,  machine  time  required  to  compute  the  isoclines 
and  trajectories  of  a  second-order  nonlinear  differential 
equation  was  about  60  sec.   This  is  an  appreciable  saving 
compared  with  hand  labor. 
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APPENDIX 
COMPUTER  PROGRAMS 
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